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ABSTRACT 


This  report  summarizes  "work  on  a  number  of  different 
but  related  topics,  as  follows:  The  effect  of  seismic  source 
depth  on  Rayleigh  wave  spectra  is  examined  for  a  dissipative  half 
space  and  for  an  elastic  layer  overlying  an  elastic  half  space. 

In  the  two-layer  earth  model,  the  higher-mode  amplitude  spectra 
Increase  relative  to  the  fundamental  as  the  source  approaches  the 
base  of  the  crust.  The  exact  ray  theory  and  matrix  method  are  in 
principle  devices  for  obtaining  the  response  of  a  layered  system 
over  a  limited  time  Interval.  Both  methods  prove  to  be  inefficient 
from  the  standpoint  of  automatic  computation  but  useful  in 
analyzing  certain  general  properties  of  the  response  functions  in 
multilayered  systems.  Geometric  ray  theory  is  used  to  study  the 
effects  of  layer  thickness  and  range  on  the  refracted  arrival 
along  a  high-speed  layer  embedded  in  an  infinite  medium.  When 
the  layer  is  thick  compared  with  the  dominant  wavelength,  the 
refracted  arrival  exhibits  range-limited  shingling  and  may  consist 
of  two  or  more  events  separated  by  equal  time  intervals  which 
depend  only  on  the  layer  thickness.  The  reflection  of  a  plane 
compressional  wave  at  a  plane  Interface  is  analyzed  with  particular 
emphasis  on  the  equation  for  continuity  of  the  Instantaneous 
energy  flux.  Inside  the  critical  angle,  Knott's  equation  gives 
continuity  of  the  instantaneous  flux,  but  beyond  the  critical 
angle  Knott's  equation  must  be  replaced  by  three  separate  conditions 


The  partition  of  energy  among  the  reflected  and  transmitted  waves 
and  the  relative  pheBes  are  tabulated  for  over  2000  cases  to  show 
the  effect  of  systematically  varying  the  comprosslonal  velocity 
ratio,  the  Poisson’s  ratios  and  the  density.  Deconvolution  tech¬ 
niques  arc  applied  to  the  body  wave  phases  generated  by  two 
explosions  and  an  earthquake  in  an  effort  to  reduce  the  distortion 
introduced  by  the  transmission  and  recording  systems. 
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iirrnoDucTiON 


This  report  summarizes  research  carried  out  at  the 
California  Research  Corporation  during  the  contract  period, 
Feoruarjf  15,  1961  through  June  l4,  1964.  The  primary  purpose  of 
this  research  was  to  investigate  theoretically  the  feasibility  of 
Improving  our  ability  to  determine  seismic  energy  source  depth. 

In  the  past  source  depth  determinations  have  been  based 
exclusively  on  the  use  of  body  waves.  However,  a  large  part  ol 
the  energy  emitted  by  a  source  in  the  crust  is  trapped  in  the 
crust  and  goes  Into  surface  waves.  Surface  v^aves  attenuate  leas 
rapidly  than  body  waves  because  they  spread  in  two  dimensions 
rather  than  in  three.  These  facts  suggest  that  at  the  stations 
closest  to  an  underground  explosion  (say  500  km),  most  of  the 
energy  may  be  in  the  form  of  surface  waves.  At  great  distances 
the  surface  waves  produced  by  a  low-yield  explosion  are  weak  com¬ 
pared  to  the  body  waves  probably  because  most  of  the  surface  wave 
energy  is  concentrated  at  higher  frequencies  which  are  scattered 
by  lateral  inhomogeneities  in  the  crust.  We  have  studied  the 
effect  of  source  depth  on  the  Rayleigh  wava  spectra  for  a  dissi¬ 
pative  half  space  and  for  an  elastic  surface  layer  overlying  an 
elastic  half  space.  Some  of  the  theoretical  predictions  have 
been  compared  with  results  from  two-dimensional  seismic  model 
studies . 
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Two  methods  were  considered  for  computing  tiid  total 
response  of  a  layered  system:  (a)  the  exact  ray  theory  and 
(b)  the  matrix  method,  m  the  exact  ray  theory  the  total  response 
la  expanded  in  an  infinite  series.  Each  term  In  the  series  des¬ 
cribes  the  response  produced  by  energy  which  reaches  the  receiver 
by  a  certain  generalized  transmission  path.  The  response  over  a 
finite  time  interval  la  determined  by  a  finite  number  of  terms  in 
the  expansion.  For  the  axially  symmetric  case,  each  term  is  given 
by  a  contour  Integral.  The  time  (and  consequently  the  cost) 
required  to  perform  the  numerical  integrations  was  found  to  be 
prohibitive  even  on  a  machine  like  the  IBM  7094. 

The  theoretical  analysis  of  the  general  term  in  the  exact 
ray  theory  expansion  revealed  that  the  response  function  for  each 
path  diverges  in  the  long-time  limit  and  that  the  rate  of  diver¬ 
gence  increases  as  the  number  of  reflections  increases.  This 
suggests  that  if  the  individual  response  functions  are  added 
together  in  an  indiscrlrplnate  way  it  is  quite  possible  for  the 
magnitude  of  the  result  at  an  Intermediate  stage  of  computation 
(before  all  the  tenris  are  computed)  to  exceed  the  magnitude  of 
the  correct  total  result  by  a  number  which  exceeds  the  number  of 
significant  figures  retained  by  the  nkiChine.  A  prescription  is 
given  for  grouping  the  ray  paths  together  in  such  a  way  that  each 
group  r..iponse  function  converges  in  the  long-time  limit. 

The  first  and  second  terms  are  derived  in  a  high- 
frequency  asymptotic  expansion  of  the  general  term  in  the  exact 
ray  theory  expansion  for  the  multi-layer  system.  The  high- 
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frequency  terms  describe  re ys  which  traverse  least-time  reflection 
and  refraction  paths.  The  factors  which  describe  the  attenuation 
produced  by  geometrical  spreading  can  be  used  in  any  system  in 
v/hich  the  radiation  field  is  axially  symmetric.  All  the  informa¬ 
tion  required  to  construct  a  high-frequency  theoretical  seismo¬ 
gram  for  a  layered  system  is  given. 

The  matrix  method  was  used  to  attack  the  multilayer 

problem  along  different  Dines.  The  matrix  method  gives  the  kernel 

» 

in  a  double  integral  transform  for  the  response  of  a  layered 
system.  The  kernel  is  a  function  of  frequency  and  phase  V'locity. 
We  restrict  our  attention  to  phase  velocities  which  exceed  the 
highest  phase  velocity  for  unattenuated  normal  mode  propagation. 

In  this  restricted  domain  the  kernel  is  non-singular.  The  trun¬ 
cated  integral  describes  in  an  approximate  v/ay  the  response  in 
the  time  interval  preceding  the  arrival  of  the  normal  modes  and 
in  this  respect  determines  a  theoretical  refraction  seismogram. 

The  feasibility  of  performing  tne  double  integration  numerically 
is  considered.  The  formal  expression  for  the  kernel  is  analyzed 
to  demonstrate  the  existence  or  absence  of  unatte.nuated  normal 
modes  in  different  models. 

.-.n  accurate  knowledge  of  the  velocity  structure  in  the 
.rust  is  required  in  the  determination  of  epicentral  position  and 
focal  depth,  in  the  use  of  equalization  techniques  to  remove  from 
surface  vjavco  the  phase  distortion  introduced  by  the  transmission 
path,  and  to  remove  from  body  v/aves  the  distortion  introduced  by 


reverberation  In  the  crust.  For  these  and  other  reasons,  the 
Branch  of  Crustal  Studies  of  the  U.  S.  Geological  Survey  has 
undertaken  an  extensive  refraction  program  to  delineate  the 
velocity  structure  in  the  crust.  One  difficulty  with  the  refrac¬ 
tion  metnod  is  that  it  cannot  detect  the  presence  of  low-velocity 
layers.  This  problem  could  be  circumvented  if  it  were  possible 
to  determine  the  base  (or  thickness)  of  the  high  velocity  zones. 
The  effect  of  layer  thickness  on  the  refracted  arrival  from  a 
layer  embedded  in  an  infinite  medium. is  investigated  in  detail. 

The  two-dimensional  seismic  model  was  used  to  study  the  refracted 
arrlvals  in.  two-  and  four-layer  models. 

The  subject  of  the  reflection  and  transmission  of  a 
plane  compressional  wave  at  a  plane  interface  is  an  old  one. 

Even  so  some  new  results  are  presented  in  this  report.  On 
physical  grounds  we  know  that  the  normal  component  of  the  instan¬ 
taneous  energy  flux  must  be  continuous  across  an  Interface. 

Inside  the  critical  angle  Knott's  equation  is  the  correct  mathe¬ 
matical  statement  of  continuity.  Outside  the  critical  angle 
Knott’s  equation  requires  continuity  of  the  net  flux  (l.e.,  the 
instantaneous  flux  integrated  over  a  period),  but  not  continuity 
of  the  instantaneous  flux.  The  mathematical  statement  of  con¬ 
tinuity  of  the  instantaneous  flux  is  derived.  The  fraction  of 
the  incident  energy  which  goes  into  the  reflected  and  tran3mltte<l 
waves  and  the  relative  phases  for  the  vertical  component  of  dis¬ 
placement  are  tabulated  for  over  two  thousand  different  combina¬ 
tions  of  the  compressional  velocity  ratio,  of  Poisson's  ratio,  and 
of  density  ratio. 
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The  earth  transmission  system  and  the  recording  system 
introduce  a  certain  amount  of  distortion  in  the  body  wave  phases 
which  lengthen  the  phases  in  time  and  reduce  resolving  power. 
Deconvolution  techniques  are  applied  to  body  wave  phases  from  two 
explosions  and  an  earthquake  in  an  effort  to  remove  the  distortion 
introduced  by  the  transmission  path. 

Each  section  of  this  report  is  complete  in  itself  and  is 
independent  of  the  preceding  and  following  sections .  The  figures 
for  each  section  are  placed  either  at  the  end  of  the  section  or 
appear  on  the  page(s)  following  fir'st  mention.  All  the  references 
are  grouped  together  at  the  end  of  the  report.  Except  for  the 
sections  on  seismic  modeling  and  inverse  filtering,  the  material 
is  covered  in  greater  detail  in  the  scientific  reports  which  have 
been  issued  during  the  course  of  the  contract. 


I.  SIGNAL-TO -NOISE  RATIO  AND  SPECTRA  OF 
EXPLOSION -GENERATED  RAYLEIGH  WAVES  IN  A 
DISSIPATIVE  HALF  SPACE 

This  preliminary  theoretical  study  is  basically  concerned 
with  the  estimation  of  the  depth  of  an  explosive  source  from  a 
knowledge  of  the  radiated  Rayleigh  wave.  The  differentiation  between 
artificial  and  natural  sources  on  the  basis  of  source  depth  requires 
a  precise  determination  of  source  depth.  Present  techniques  for 
determining  source  depth  at  short  ranges  rely  exclusively  on  the 
use  of  body  wave  phases.  These  methods  presuppose  good  control 
at  short  ranges.  This  may  be  difficult  to  obtain.  It  is  reasonable 
to  inquire  whether  accurate  source  depth  information  can  be 
extracted  from  surface  waves  at  moderate  ranges,  say,  of  the  order 
of  500  kms.  At  great  distances  the  surface  waves  generated  by 
low-yield  nuclear  explosions  are  generally  not  detected  probably 
because  most  of  the  energy  goes  into  the  higher  frequency  components 
which  are  scattered  by  lateral  Inhomogeneltles  in  the  upper  crust. 

To  investigate  the  effect  of  source  depth  on  the  surface 
wave  spectrum  and  the  factors  which  govern  our  ability  to  recover 
that  spectrum,  we  consider  the  Rayleigh  wave  generated  by  an 
explosion  in  a  homogeneous  and  Isotopic  half  space.  Thjs  is  the 
simplest  medium  which  permits  the  propagation  of  Rayleigh  waves. 

The  explosive  charge  is  assumed  to  be  situated  at  the  center  of  a 
spherical  cavity  of  radius  a.  In  what  follows,  decoupling 
(Latter,  Le  Levier,  Martinelli,  and  McMillan,  I96I)  is  assumed 
throughout . 
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The  Laplace  transform  of  the  radial  stress  in  a 
spherically  symmetric  compressional  wave  is 


^rr=-^ 


_sr/oc_ 


V  oc^vsry  ac^W  J  > 


where  T  is  the  radial  distance  from  the  center  of  the  shot  cavity, 
S  is  the  Laplace  transform  variable,  0^  is  the  density,  /6  is  the 
shear  velocity,  and  OL  is  the  compressional  velocity.  Following 
Latter,  we  approximate  the  stress  at  the  cavity  wall  by  a  step. 
This  requirement  determines  in  the  form: 


where  P  is  the  magnitude  of  the  step  in  stress.  Because  of 
preferential  absorption  of  the  higher  frequencies,  only  the  wave¬ 
length  components  which  are  large  compared  to  the  cavity  radius  are 
important  at  distance.  The  assumption  |c(/Scl|»  1  ,  reduces  (2) 

(3) 

To  determine  whether  (3)  is  valid,  we  must  compare  the  shortest 
wavelengths  in  the  Rayleigh  wave  train  at  500  km  with  the  minimum 
cavity  radius  required  to  achieve  decoupling.  The  quantity  PCL^ 
determines  the  strength  of  the  source  and  is  related  to  the 
yield,  W,  by 

where  'Y  is  the  ratio  of  the  specific  heats  of  the  gases  in  the 
cavity. 
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Sherwood  and  Spencer  (1962)  have  derived  an  expression 
for  the  spectrum  of  the  vertical  component  of  the  particle  velocity 
in  the  Rayleigh  wave  at  the  free  surface.  The  phase  spectrum  is 


(5) 


where  -f  is  the  frequency,  p  is  the  range,  and  C  is  the  Rayleigh 
wave  velocity.  It  is  to  be  noted  that  the  phase  spectrum  is 
independent  of  the  source  depth,  H.  The  amplitude  spectrum  is 


°  pi/s 


3^ 


(6) 


where  M  is  the  TNT  mass  equivalent  in  kilotons  and 

TC  _Ca6)lo'°Bc'^^ 

B  is  a  complicated  function  of  Poisson’s  ratio  which  decreases 
monotonically  from  about  1.2  to  0.2  as  the  Poisson’s  ratio 
increases  from  zero  to  one -half.  It  is  evident  from  (6)  that 
the  higher  frequency  components  decay  more  rapidly  with  source 
depth  than  the  lower  frequencies.  Hence,  the  amplitude  spectrum 
of  the  Rayleigh  wave  becomes  relatively  richer  in  high  frequencies 
as  the  source  depth  H  is  diminished. 

Unfortunately,  the  absorptive  properties  of  the  medium 
have  an  effect  similar  to  source  depth  in  diminishing  the  relative 
content  of  the  higher  frequency  components.  It  seems  that  the 
loss  properties  of  a  solid  medium  can  often  be  approximately  repre¬ 
sented  by  a  mechanical  quality  factor,  Q,  which  does  not  vary 
much  with  frequency  (see  O’Brien  (I96I)  for  bibliography). 
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Hueter  and  Bolt  (1955)  state  that  the  absorption  of  pro¬ 
gressive  waves  Is  approximately  expressible  as 


T-,  ^ 


(7) 


An  interpolation  of  measurements  in  the  vicinity  of  .01  to  30  cps 
(see  O’Brien  (1961)  and  Press  (195^))  indicates  that  in  the 
frequency  range  of  Interest  (l  to  .1  cps),  the  Q  for  Rayleigh 
waves  should  be  very  approximately  100.  It  does  seem,  however, 
that  lower  frequency  Rayleigh  waves  have  a  higher  Q  value.  Also, 
the  higher  frequencies  generally  experience  considerable  scatter¬ 
ing  In  the  Inhomogeneous  near-surface  region  which  should  result 
in  a  decrease  in  the  effective  Q  value.  Hence,  the  effective 
Q  for  a  surface  wave  mode  should  tend  to  decrease  with  increasing 
frequency,  but  we  are  not  presently  in  a  position  to  give  precise 
quantitative  values.  This  is  important,  for  we  see  from  (6)  and 
(7)  that  the  source  depth  and  dissipation  have  a  similar  frequency- 
dependent  effect  on  the  amplitude  spectrum  of  the  Rayleigh  wave. 

It  is  thereby  apparent  that  an  error  in  our  knowledge  of  Q  will 
limit  our  ability  to  estimate  the  source  depth  accurately. 

The  combined  effects  of  source  depth  and  dissipation 
determine  the  amplitude  spectrum  in  the  form 


Aft')- -  DLM 

pl/2 

where  • 

Assuming  that  the  velocities  C  and  cL  are  known  accurately,  we 


(8) 


(9) 


see  that  the  incremental  changes  in  the  other  parameters  are 
related  by 
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v/here  {ix )  represents  a  small  Increment  in  X  •  Even  under  the 
assumption  that  ^  can  be  accurately  estimated  from  the  ampli¬ 
tude  spectrum  of  uhe  Rayleigh  wave  train,  and  even  if  the  range 
p  is  known  accurately,  it  is  evident  that  there  is  an  error  in 
source  depth,  given  by 


(10) 


To  appreciate  the  magnitude  of  this  error,  let  us  enter  plausible 
values  into  (lO).  For  a  resonable  value  of  C/oit  a  range  p  of 
500km#  and  an  error  (  <$”  Q)  of  10  in  a  Q  of  about  100,  we  have 

0.3  Km. 

This  is  about  the  resolution  in  '•"'urce  depth  that  we  are  seeking, 
but  it  seems  that  it  will  take  considerable  effort  to  attain  it. 
The  normalized  Rayleigh  spectra  defined  by 

uc-f)=  p'*/l(-f)/M 

are  plotted  in  Figure  1. 

Our  ability  to  extract  information  from  the  Rayleigh 
wave  signal  is  limited  by  the  microseism  or  noise  background. 

We  assume  that,  at  the  time  and  place  of  detection,  the  charac¬ 
teristics  of  the  background  noise  N(.p)  are  similar  to  those  of 
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the  minimum  noise  curve  depicted  by  Brune  and  Oliver  (1959)  and 
reproduced  in  Figure  1.  Accordingly,  our  final  estimates  refer 
to  approximately  the  maximum  signal -to -noise  ratio  that  we  can 
expect  for  a  decoupled  shot. 

Actual  seismograms  are  the  result  of  subjecting  the 
true  ground  motion  to  some  composite  filter  system.  This  filter 
system  is  usually  a  reasonably  simple  combination  of  a  conven¬ 
tional  seismometer  and  galvanometer.  However,  it  could  in 
principle  be  a  far  more  sophisticated  filter  which  one  might 
obtain  by  means  of  analog  or  digital  operations  on  the  output 
of  a  seismometer.  The  Importance  of  recording  a  wide  range  of 
seismic  frequencies  with  good  accuracy  in  conjunction  with  ade¬ 
quate  filtering  operations  cannot  be  overemphasized.  Sherwood 
and  Spencer  (1962)  consider  two  filter  systems:  (l)  the  Benioff 
short  period  system  installed  at  Wichita*  and  (2)  a  more  nearly 
optimum  system  in  which  the  filter  characteristic  is  a  function 
of  the  spectral  content  of  both  the  Rayleigh  signal  and  the 
ambient  ground  motion. 

As  our  definition  of  signal -to-nolse  ratio,  we  adopt 
the  following  form: 


♦The  transmission  curve  for  particle  velocity  for  the  Wichita 
Instrument  is  given  in  the  article  by  Gudzin  and  Hamilton  (196I), 
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Average  level  of  signal  amplitude  spectrum _ 

Average  level  of  conflicting  noise  amplitude  spectrum 

Recording  the  signal  and  noise  through  the  short  period  Benioff 
system  gives  a  signal -to -noise  ratio  per  kiloton  of  TNT  |ls/N)/^ 
which  varies  with  range  as  shown  in  Figure  2. 

It  must  be  emphasized  that  our  results  are  valid  only 
for  the  decoupled  source,  step  function  in  pressure  input,  con¬ 
stant  Q  dissipation,  half-space  model,  and  noise  structure  which 
we  have  assumed.  A  coupled  shot  is  far  more  efficient  in  gen¬ 
erating  seis''*G  waves  than  a  decoupled  shot,  and  its  spectrum 
may  be  different.  It  is  particularly  important  to  compare  the 
long  wave  length  content  of  the  two  spectra  in  situations  where 
the  free  surface  or  some  other  major  impedance  discontinuity 
lies  within  the  nonlinear  zone.  In  a  layered  earth  the  Rayleigh 
waves  propagate  in  more  than  one  mode,  and  each  mode  exhibits 
dispersion.  Dispersion  causes  the  Rayleigh  wave  to  spread  out 
over  a  relatively  long  time  interval  with  a  subsequent  decrease 
in  the  signal -co-noise  ratio.  When  the  phase  velocity  vari¬ 
ation  with  freauency  is  known,  it  is  possible  partially  to 
compensate  for  the  deleterious  effect  of  dispersion  by  a  tech¬ 
nique  such  as  that  due  to  Akl  (i960).  The  variation  of  absorp¬ 
tion  with  depth  and  the  lateral  inhomogeneity  of  the  wave  guide 
probably  cause  the  Q  to  decrease  with  increasing  frequency. 

These  factors  must  be  considered  in  comparing  our  results  with 
observed  signal-to-noise  ratios. 


A(f)  IS  IN  METERS i  P,  METERS  i  M,  KILOTONSJ 


FIGURE  I 

NORMALIZED  RAYLEIGH  WAVE  AND 
NOISE  AMPLITUDE  SPECTRA 


LE  19-293 


NOISE  VELOCITT  (METERS/SEC.)  N(f) 


FIGURE  2 

APPROXIMATE  SIGNAL  TO  NOISE  RATIO 
PER  KILOTON  OF  TNT  AS  A  FUNGI lON 
OF  RANGE,  DETECTOR  IS  A  SHORT 
PERIOD  RENIOFF  INSTRImENT.  H  =  0. 


LE  19-296 
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li.  THE  EFFECT  OF  SOURCE  DEPTH  ON  RAYLEIGH  WAVE  MOTION 

IN  A  LAYERED  EARTH 

In  comparison  with  the  more  conventional  use  of  the 
body-wave  phases,  the  use  of  surface  waves  for  source  depth 
determination  has  both  advantages  and  disadvantages .  A  prominent 
disadvantage  of  surface  waves  for  source  depth  woi'k  is  their 
longer  periods.  It  is  evident  that  the  effect  of  small  changes 
in  source  depth  will  be  reflected  only  in  the  higher  frequency 
portions  of  the  wave  train.  The  results  of  our  Rayleigh  wave 
work  substantiate  this  point.  For  the  case  of  an  earth  model 
consisting  of  a  30  kilometer  crust  overlying  a  half-space  of  mantle 
material,  t..e  effect  of  a  change  in  depth  for  a  source  within 
the  crust  is  apparent  only  at  periods  less  than  10  seconds. 

A  possible  advantage  of  surface  waves  in  seismic  studies 
is  that  they  undergo  less  attenuation  due  to  geometrical  spread¬ 
ing  than  body  waves.  In  addition,  surface  waves  are  less  affected 
by  small  changes  in  velocity  structure  than  body -wave  phases. 

The  use  of  body-wave  arrivals  for  source  depth  determination 
requires  an  accurate  knovjle.'ge  of  the  velocity  structure  in  the 
region  between  the  source  and  the  receiver.  The  recognition  of 
source  depth  effects  on  surface  'wave  motion  does  not  require  as 
detailed  a  knowledge  of  the  velocity  structure. 
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The  general  method  of  solution  was  as  follows: 
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h 
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5  source 
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^i>A>  pi 


.recelYer 


^2  ^^2)  p2. 

A  point  impulsive  source  is  located  at  a  depth, h  ,  within  a 
surface  layer  of  thickness  H  overlying  a  semi-infinite  half -space. 
The  motion  is  recorded  at  a  horizontal  offset /L  from  the  source. 
cLi  ^  )  and  p  are,  respectively,  compressional  wave  speed,  shear 
wave  speed,  and  density.  The  problem  is  to  obtain  an  expression, 
capable  of  being  evaluated,  for  the  motion  at  a  receiver  at 
distance  A,.  It  develops  that  we  cannot  obtain  a  closed  algebraic 
expression  for  the  total  motion  recorded  from  the  point  source. 

The  best  we  can  do  is  to  obtain  an  integral  solution  for  the 
motion.  It  is  of  the  form 
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It  is  not  possible  to  evaluate  this  integral  exactly,  but  various 
approaches  have  been  used  to  obtain  certain  parts  of  the  solution. 
One  approach  (Pekeris,  19^8;  Spencer,  i960)  is  to  expand  the 
integrand  in  an  infinite  series,  each  term  of  which  represents 


the  contribution  from  energy  that  has  traveled  a  particular  path 
to  the  receiver.  Each  of  these  new  integrals  can  be  evaluated  by 
Cagniard's  method.  A  second  approach,  and  the  one  used  in  the 
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Rayleigh  work,  is  to  attempt  a  direct  evaluation  of  the  integral 
by  means  of  contour  integration.  In  such  an  approach  the 
Rayleigh  wave,  or  normal  mode  motion,  arises  from  the  contri¬ 
butions  of  the  residues  at  the  real  poles  of  the  integrand. 

The  contributions  of  these  residues  were  evaluated  numerically 
for  the  present  work.  The  solutions  obtained  gave  the  spectrum 
of  the  normal  mode  motion  for  both  the  horizontal  and  vertical 
components  of  displacement  resulting  from  a  buried  point  source. 

For  the  case  in  which  the  model  consists  of  a  half¬ 
space  of  uniform  material,  the  effect  of  placing  the  source 
deeper  beneath  the  free  surface  is  well  known  (Rayleigh  1885). 

The  effect  of  source  depth, Z  ,  upon  the  spectrum  is  given  by  a 

c  1 

term  0  ,  where  t  is  frequency  and  0  is  a  term  independent  of 

source  depth.  From  this  it  is  evident  that  a  deeper  sou"*ce 
results  in  longer  period  motion.  For  the  case  of  a  source  within 
a  layered  elastic  medium,  the  effect  of  source  depth  upon  the 
normal  mode  motion  is  not  so  simple.  In  such  a  system,  an  infinite 
number  of  different  modes  of  propagation  are  theoretically 
possible,  and  the  total  observed  motion  is  the  sum  of  an  Infinite 
number  of  these  modes.  The  frequency  spectra  and  the  relative 
excitation  of  the  modes  are  functions  of  source  depth.  Hence  the 
total  observed  motion  will  also  be  a  function  of  source  depth. 

The  aim  of  our  work  on  normal  mode  motion  was  to  determine 
explicitly  just  how  the  source  depth  does  affect  the  total  motion, 
at  least  for  the  first  few  modes. 
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The  first  case  considered  was  that  of  a  liquid  layer 
overlying  a  liquid  half -space.  This  case,  which  is  algebraically 
much  simpler  than  the  solid  case  considered  later,  served  as  an 
illustration  of  the  method  and,  in  addition,  gave  a  partial 
check  on  the  answers,  since  the  case  had  been  previously  consid¬ 
ered  by  Pekeris  (19^8)  and  certain  of  his  solutions  can  be  applied. 
The  expression  for  the  horizontal  and  vertical  components  of 


displacement,  (j^  and  lY  respectively,  are  given  by 


2  JjCKA)clK  ^ 


W — 2  K^(xi«.T'Z)3;(K;i)dK. 


A  is  the  Rayleigh  determinant,  K  is  the  wave  number, 

_  '  oLf  / 

and  “V^i  =  K€|  .  is  a  3  by  3  determinant  involving  various 

parameters  of  the  system.  Evaluation  of  these  Integrals  in  the 


complex  zeta  plane,  r-K+iT  ,  yields  the  following  expressions 
for  the  normal  mode  portions  of  the  solution: 


W— 2#if  X 


represents  dA/ d  K 


The  summation  is  over  the  real  poles 


of  the  Integrand. 
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The  expressions  for  and  W  were  evaluated  numerically 
on  an  IBM  7090  digital  computer.  The  programs  were  written  for 
variable  values  of/l,Z,h  ,  H,  Ct  ^  /3  ,  and  p  .  Figures  1  through 
6  show  examples  for  the  effect  of  source  depth  on  the  vertical 
component  of  particle  displacement.  Figure  1  is  a  plot  of  the 
individual  amplitude  spectra  for  the  first  six  normal  modes  for 
a  source  depth  of  1/6  the  layer  thickness.  Figure  2  shows  the 
amplitude  spectra  for  a  source  depth  of  29/30  the  layer  thickness. 
Figures  3  through  6  show  the  amplitude  spectra  for  the  motion 
consisting  of  the  sum  of  the  first  six  normal  modes  as  a  function 
of  source  depth.  The  receiver  offset  is  l8.3  times  the  layer 
thickness,  and  four  values  of  source  depth  are  used. 

The  general  conclusions  for  the  liquid  case  are  that 
the  higher  modes  are  important  and  give  significant  contributions 
to  normal  mode  motion  in  a  liquid  system.  The  amplitude  spectra 
of  the  individual  normal  modes  are  broad  and  contribute  motion 
over  a  large  range  of  H  above  the  cut-off  frequency.  The  spectra 
for  the  total  summed  motion  in  the  liquid-liquid  case  are  quite 
complicated  because  of  the  contribution  of  the  many  modes  that 
are  required  to  describe  the  total  motion.  The  effect  of  source 
depth  on  the  finer  structure  of  the  spectra  for  the  total  normal 
mode  motion  in  the  liquid-liquid  case  most  probably  will  not  be 
observed  because  of  the  complicated  nature  of  the  summed  spectra. 
The  most  important  feature  seems  bo  be  a  broadening  of  the  spec¬ 
trum  as  the  source  is  placed  deeper  within  the  layer.  Longer 
period  components  are  relatively  stronger  for  deeper  sources. 
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For  the  case  of  a  solid  elastic  layer  overlying  a  solid 
half -space,  equations  similar  to  those  for  the  lic.uid  case  were 
obtained.  The  solid  case  equations  are,  however,  quite  lengthy 
and  v;ill  not  be  written  down  here.  It  was  found  that  the  higher 
normal  modes  are  not  nearly  so  important  in  the  solid  case  as 
they  are  for  the  liquid  case.  Three  modes  are  sufficient  to  des¬ 
cribe  the  total  Rayleigh  wave  motion  with  reasonable  accuracy. 

The  spectra  of  the  Individual  normal  modes  are  more  sharply  peaked 
than  in  the  liquid  case.  Figures  7  through  9  show  examples  of 
the  amplitude  spectra  for  the  horizontal  component  of  displacement 
of  the  first  three  Rayleigh  modes.  Examples  are  taken  for  source 
depths  of  1/30  the  layer  thickness,  1/3  the  layer  v._.ickness,  and 
2/3  the  layer  thickness.  Figures  10  through  15  show  the  spectrum 
of  the  total  motion  resulting  from  the  sum  of  the  first  three 
normal  modes.  This  again  is  for  the  horizontal  component  of  dis¬ 
placement,  and  the  source  depths  are  taken  from  near  the  free  sur¬ 
face  at  1/30  the  layer  thickness  to  near  the  layer-half  space 
interface  at  29/30  the  layer  thickness.  The  effect  of  source 
depth  on  the  spectrum  is  apparent  only  where  the  contribution  of 
the  second  mode  begins  to  be  significant.  This  is  at  a  value  of 
■fH/o(.|  =  .5  for  a  30-kilometer  crust,  and  for  this  case  corresponds 
to  a  period  of  10  seconds.  An  appreciable  difference  in  the  spectra 
as  a  function  of  source  depth  appears  only  for  periods  of  7  seconds 
or  less.  The  peak  of  the  Rayleigh  wave  motion  for  a  source  within 
a  30-kilometer  crust  is  at  a  period  of  20  seconds.  A  source  near 
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the  free  surface  generates  predominantly  first  mode  motion.  As 
the  source  Is  placed  deeper  within  the  layer,  the  higher  modes, 
and  higher  frequencies,  become  relatively  more  important. 

In  summary,  there  are  theoretically  predicted  differ¬ 
ences  in  the  amplitude  spectrum  of  Rayleigh  wave  motion  as  a 
function  of  source  depth.  These  differences  are  evident  only 
in  the  higher  frequency  portion  of  the  wave  train  at  periods 
less  than  10  seconds.  V/hether  these  effects  will  be  useful  in 
practice  in  the  pi esence  of  microseisra  noise  is  a  question  best 
ansv/ered  by  an  analysis  of  field  seismograms. 
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FIGURE  4 

SPECTRUM  OF  SUM  OF  FIRST 

SIX  NORMAL  MODES. 
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FIGURE  6 


SPECTRUM  OF  SUM  OF  FIRST 

SIX  NORMAL  MODES. 
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FIGURE  7 

INDIVIDUAL  SPECTRA  FOR  HORIZONTAL  COMPONENT  OF 
DISPLACEMENT  FOR  FIRST  THREE  RAYLEIGH  MODES. 
SOLID  LAYER  OVER  SOLID  HALF  SPACE. 
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FIGURE  8 
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INDIVIDUAL  SPECTRA  FOR  HORIZONTAL  COMPONENT  OF 
DISPLACEMENT  FOR  FIRST  THREE  RAYLEIGH  MODES. 
SOLID  LAYER  OVER  SOLID  HALF  SPACE. 
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FIGURE  9 

INDIVIDUAL  SPECTRA  FOR  HORIZONAL  COMPONENT  OF 
DISPLACEMENT  FOR  FIRST  THREE  RAYLEIGH  MODES, 
SOLID  LAYER  OVER  SOLID  HALF  SPACE. 
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SPECTRUM  OF  SUM  OF  HORIZONTAL  OlSi  ^ACEMENTS 
FOR  FIRST  THREE  RAYLEIGH  MODES. 
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FIGURE  II 


SPECTRUM  OF  SUM  OF  HORIZONTAL  DISPLACEMENTS 
FOR  FIRST  THREE  RAYLEIGH  MODES. 
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FIGURE  13 
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SPECTRUM  OF  SUM  OF  HORIZONTAL  DISPLACEMENTS 

FOR  FIRST  THREE  RAYLEIGH  MODES. 
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III.  SEISMIC  MODELING 


A.  Explosive  Sources 

Initially  the  model  work  was  set  up  to  supply  corrobora¬ 
ting  data  for  our  theoretical  work  dealing  with  the  effect  of 
source  depth  upon  Rayleigh  wave  motion.  The  theoretical  solutions 
give  the  motion  produced  by  an  impulsive,  point,  compressional 
source.  Because  of  this  we  wanted  a  source  with  as  wide  a  band 
as  possible  and  hence  began  our  work  using  an  explosive  source. 

In  this  section  we  review  the  results  which  were  obtained  using 
this  source. 

Two-dimensional  seismic  models  were  used  for  all  the 
cases  which  we  studied.  It  is  easier  to  construct  and  work  with 
two-dimensional  models  and  then  extend  the  results  to  three- 
dimensional  wave  propagation  than  it  is  to  construct  three- 
dimensional  models.  The  basic  model  consisted  of  a  sheet  of 
cold-rolled  steel,  8  feet  long,  4  feet  wide,  and  .060  inch  thick. 

For  the  wavelengths  used  in  the  m.odeiing  experiments,  the  thick¬ 
ness  of  this  sheet  is  vanishingly  small.  A  layer  having  elastic 
properties  different  from  those  of  the  steel  was  formed  along 
one  edge  of  the  sheet  by  first  milling  out  a  groove  .030  inch 
deep,  1  inch  wide,  and  8  feet  long.  This  mllled-out  section  was 
then  filled  with  a  low  velocity  material  to  simulate  the  case  of 
a  low-velocity  layer  overlying  a  high-velocity  half-space.  In 
the  earlier  work  the  groove  was  filled  with  a  ^0-^0  lead-tin  solder. 
Unfortunately,  it  was  difficult  to  hold  the  thickness  of  the  solder 


layer  to  a  proper  tolerance  without  a  further  milling  operation. 

A  different  method  was  employed  for  latei  models.  This  method 
consisted  of  cementing  a  copper  strip  .030  inch  thick,  1  inch 
wide,  and  8  feet  long  into  the  milled -out  section.  This  gives 
good  control  over  the  thickness  of  the  composite  layer  hut  intro¬ 
duces  the  problem  of  obtaining  good  coupling  betv/een  the  copper 
and  the  steel.  A  semi-rigid  epoxy  cement,  Eccobond  45,  was  used 
in  constructing  the  models  for  the  explosive  source  work.  A 
slightly  flexible  epoxy  was  chosen  lo  triat  the  bond  between  the 
copper  and  the  steel  might  withstand  a  small  amount  of  flexing 
^'fithout  breaking  and  so  that  the  model  would  be  better  able  to 
resist  the  high  local  stresses  resultlr.g  from  the  explosions. 

Th:'3  type  of  epoxy  did  give  a  somewhat  smaller  signal  amplitude 
than  a  similar  model  .vitii  a  rigid  epoxy,  but  this  was  not  a 
problem  sine  ■  there  was  always  sufficient  signal  strength  from 
Che  explosive  source. 

Figure  1  shows  the  construction  of  the  twe  -dimensional 
seismic  models.  Wave  propagation  in  two-dimensional  seismic  models 
has  been  discussed  by  Oliver,  et  al.,  (195-0-  Compresslonal  v/ave 
energy  propagates  with  the  plate  wave  velocity  of  the  material. 

The  shear  wave  energy  propagates  with  the  s lear  velocity  in  an 
-infinite  solid  and  is  unaffected  by  the  two-dimensional  nature  of 
the  system.  Paylelgh  waves  in  two  dimensions  propagate  in  a 
slightly  different  manner  from  those  in  the  infinite  solid  case. 

The  major  difference  between  the  infinite  solid  and  the  two- 
dimensional  case  is  that  the  plate  dilate tional  velocity  replaces 
the  bulk  P  wave  velocity  in  the  Rayleigh  equation. 


CONSTRUCTION  OF  THE  2  DIMENSIONAL  SEISMIC  MODELS 


The  seismic  velocities  in  the  composite  layer  are 
governed  by  the  relative  thicknesses^  densities,  and  elastic  para¬ 
meters  of  the  materials  making  up  the  laytj.  For  wavelengths  which 
are  long  compared  to  the  total  thickness  of  the  layer,  the  body- 
wave  velocities  in  the  composite  system  are  given  by  Angona  (i960) 
as : 

_  0iiP,CiScl;P^Ca^+  .  .  ■  ■  0i„PnC„^ 

<<|P,  +a2p2+  •  •  ■  ctnPn 

where 

Otp  is  the  fraction  of  material  n > 

p  is  the  density  of  material  n  > 

J  n 

is  the  acoustic  wave  velocity  of  material  fl  . 

The  phase  and  group  velocities  for  the  first  three 
Rayleigh  modes  for  the  model  shown  in  Figure  1  are  given  in 
Figure  2.  The  phase  velocities  of  the  higher  modes  lie  between 
the  shear  velocity  in  the  steel  and  that  in  the  composite  layer. 

The  first  mode  phase  velocities  lie  between  the  Rayleigh  velocities 
in  the  steel  plate  and  in  the  composite  layer. 

The  source  used  in  these  model  experiments  was  a  small 
charge  of  silver  acetylide.  Silver  acetyllde  is  formed  as  a  pre¬ 
cipitate  when  acetylene  gas  is  bubbled  through  an  aqueous  solution 
of  silver  nitrate.  When  dry  this  precipitate  is  a  powerful  explo¬ 
sive  which  can  be  detonated  by  either  heat  or  shock,  A  certain 
amount  of  care  is  required  in  the  preparation  and  use  cf  this 
explosive . 
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FIGURE  2 


(H  is  layer  thickness) 

PHASE  AND  GROUP  VELOCITY  FOR  FIRST  THREE  RAYLEIGH  MODES  IN  2  DIMENSIONAL 
SEISMIC  MODEL. 


Chemically  pure  silver  nitrate,  distilled  water,  and 
purified  acetylene  (99-5  per  cent  minimum  purity)  are  used  In 
making  the  explosive.  The  use  of  tap  water  can  lead  to  precipi¬ 
tates  which  are  too  unstable  to  be  handled  safely  in  model  work. 
Four  grams  of  silver  nitrate  are  dissolved  in  I50  milliliters  of 
distilled  water,  and  acetylene  is  bubbled  through  the  mixture 
until  the  yellow  precipitate  is  formed  (about  five  minutes). 

The  precipitate  is  collected  on  a  filter  paper  and  dried  by 
blotting  with  dry  filter  papers  until  it  becomes  a  plastic  mass. 
The  silver  acetyllde  charges  are  then  shaped  and  mounted  on  a 
small  rod.  Figure  3  shows  the  details  of  the  die  in  which  the 
charges  are  made.  Figure  4  shows  the  procedure  for  making  the 
charges.  The  charges  are  .048  inch  in  diameter,  0.1  inch  in 
length,  and  are  mounted  on  the  end  of  a  steel  rod.  The  charges 
are  detonated  by  means  of  an  electrically  heated  wire  which  is 
positioned  close  to  the  explosive.  The  flash  of  the  charge  is 
used  to  trigger  the  o.^  lloscope  on  which  the  model  records  are 
displayed. 

The  receiver  which  was  first  used  in  this  particular 
model  work  is  essentially  a  condenser  microphone.  The  surface 
of  the  model,  or  one  face  of  a  small  square  hole  cut  into  the 
model,  forms  one  plate  of  the  condenser  microphone.  The  other 
plate  is  mounted  on  the  end  of  a  shielded  lead  connecting  directly 
to  the  grid  of  a  cathode  follower.  From  the  cathode  follower  the 
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1.  Charge  die  by  tapping  several  times  Into  well  filled  with 
damp  explosive. 

2.  Turn  die  upside  down  over  clean  hole  and  force  excess  explosive 
out.  If  more  compaction  is  required,  this  may  flret  be  done  on 
a  smooth  surface  and  then  forced  out  into  hole. 

2A.  Cut  explosive  off  even  with  end  of  die  head. 

3.  Remove  punch,  Insert  mounting  rod,  place  die  head  on  smooth 
surface  and  tap  mounting  rod  to  seat  explosive. 

Reverse  die,  tap  end  of  mounting  rod  and  remove. 


FIGURE  A 
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output  goes  directly  to  an  oscilloscope.  Figure  5  shows  both  the 
circuit  for  the  condenser  microphone  receiver  and  the  details  of 
the  phototube  circuit  which  is  used  to  trigger  the  sweep  of  the 
oscilloscope. 

The  use  of  the  explosive  source  with  the  condenser 
microphone  receiver  gives  wide  band  records  with  clearly  defined 
arrivals.  Figure  6  shows  four  of  these  records.  The  source  is 
positioned  either  at  the  surface  or  within  a  round  hole  .064  inch 
in  diameter.  The  receiver  is  located  either  on  the  surface  or  in 
a  small  square  hole.  The  arrivr.ls  of  all  of  the  various  wave 
types  at  the  receiver  ere  clearly  defined  and  quite  sharp.  In 
spite  of  their  apparent  good  quality,  however,  these  records 
were  not  appropriate  for  our  modeling  of  the  effect  of  source 
depth  upon  model  seismograms  because  of  the  following  considera¬ 
tions.  We  were  Interested  in  modeling  a  mathematical  solution 
that  postulated  a  point  compressicnal  source  at  depth.  In  order 
to  simulate  this  source,  a  small  hole  was  drilled  in  the  model, 
and  an  explosive  charge  was  detonated  in  this  hole.  No  matter 
how  carefully  the  explosive  was  centered  in  the  hole,  a  large 
shear  wave  was  invariably  generated.  This  made  the  source 
inapplicable  for  our  work.  An  alternative  approach  to  determine 
the  effect  of  source  depth  upon  the  seismic  wave  form  received 
at  the  surface  is  to  make  use  of  a  reciprocity  relationship  between 
the  dllataticn  and  th*  vertical  component  of  displacement 
(Rayleigh,  1945).  For  a  given  separation  between  source  and 
receiver,  the  vertical  displacement  recorded  at  the  surface  from 
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a  dilatational  source  at  a  depth,  h  >  Is  the  sans  as  the  dilatation 
at  the  depth,  h  t  resulting  from  a  vertical  force  applied  at  the 

4  » 

surface.  Thus  for  the  modeling  work  we  could  apply  a  vertical 
force  (by  means  of  an  explosive  charge)  at  the  surface  of  the 
model  and  record  the  dilatation  at  depth.  This  would  be  equiva¬ 
lent  to  the  effect  of  a  change  in  depth  of  a  dilatational  (i.e,, 
pure  compresslonal)  source  on  the  vertical  component  of  displace¬ 
ment  at  the  surface.  The  method  which  was  first  employed  to  take 
advantage  of  the  above-mentioned  recrlpoclty  was  to  place  an 
explosive  source  at  the  surface  of  the  model  and  to  locate  a  con¬ 
denser  microphone  in  a  small  hole  at  a  specified  depth  in  the 
model.  Unfortunately,  this  method  failed  to  give  results  useful 
in  checking  our  mathematical  solutions  because  of  the  difficulty 
in  relating  in  a  precise  way  the  output  of  the  condenser  micro¬ 
phone  to  the  dilatation.  There  was  obviously  some  non-dllatatlonal 
deformation  of  the  hole  as  the  seismic  energy  arrived  and  this 
was  also  recorded  by  the  condenser  microphone. 

The  approach  which  was  finally  adopted  was  to  use  strain 
gages  to  record  the  dilatation  at  depth  in  the  model.  The 
receivers  were  heavily-doped  silicon  semi-conductor  strain  gages, 
type  POl -05-120  made  by  Micro  Systems,  Inc.  These  gages  are 
.060  inches  long,  .005  Inches  wide,  and  .002  Inches  thick.  At 
each  desired  depth  location  two  gages,  crossed  at  right  angles 
to  one  another,  were  cemented  to  the  model.  One  gage  then  gave 
the  JU/JjC  component  of  the  dilatation  while  the  other  gave  the 
dV/JiJ  component.  The  outputs  of  the  two  strain  gages  were 
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summed  in  series  to  give  the  total  dilatation,  0— -4-JV/oj^  • 

The  original  configuration  had  four  gages  mounted  at  each  receiver 
location,  two  on  each  side  of  the  model  exactly  opposite  one 
another.  The  purpose  of  this  arrangement  was  to  cancel  out  flexural 
waves  in  the  model.  It  developed,  however,  that  the  flexural  waves 
were  small  and  did  not  Interfere  with  the  normal  mode  motion,  so 
the  additional  two  gages  were  eliminated.  Figure  7  shows  the 
placement  of  the  gages  at  four  different  depths  on  the  model. 

The  source  is  located  as  shown  at  the  top  of  the  model.  This 
particular  position  of  the  source  largely  eliminated  flexural 
vibrations  and  gave  quite  reproducible  results.  Figure  8  gives 
an  example  of  two  shots  recorded  at  the  same  location.  The  shots 
were  at  the  surface  of  the  model  and  the  receivers  were  located 
near  the  layer-half  space  Interface  at  an  offset  of  40  centimeters. 
As  can  be  seen  from  the  example,  there  are  essentially  no 
differences  between  the  two  records. 

Figures  9  and  10  show  the  records  and  the  corresponding 
amplitude  spectra  from  four  different  receiver  locations.  The 
receivers  were  located  as  near  as  possible  to  the  top  of  the  model, 
at  a  depth  of  1/6  the  layer  thickness,  at  a  depth  2/3  the  layer 
thickness,  and  as  near  as  possible  to  the  layer-half  space  boundary. 
The  one-inch  layer  on  the  model  was  assumed  to  represent  the 
crust  of  the  earth  and  the  modeling  results  were  adjusted  to  a 
one-inch-equals -30-kllometer  change  in  scale.  The  abscissa  scale 
for  the  period  of  the  spectra  is  in  seconds.  The  amplitude  spectra 
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can  then  be  compared  with  the  results  derived  from  the  theoretical 
normal  mode  solutions.  The  modeling  results  show  that  from  a 
period  of  7  or  8  seconds  out  to  longer  period  motion,  the  spectrum 
for  the  surface  source  is  flatter  than  the  spectra  for  the  deeper 
sources.  As  the  source  is  placed  deeper  within  the  layer,  the 
spectra  become  slightly  more  Irregular  and  the  high  frequency 
components  increase  slightly  relative  to  the  longer  period  components. 
This  same  effect  was  apparent  in  the  theoretical  solutions. 

B.  Piezoelectric  Sources 

The  use  of  piezoelectric  sources  and  receivers  was  motiva¬ 
ted  by  two  basic  considerations.  In  the  first  place,  the  use  of 
explosive  sources  limits  the  choice  of  materials  from  which  the 
models  can  be  constructed.  Only  steel  was  able  to  withstand  the 
force  of  the  explosion  without  deforming,  and  even  it  would  deform 
under  repeated  shots.  Copper  and  lead  deform  easily  and  a  single 
shot  will  cause  cracking  and  spalling  in  a  Plexiglas  model.  The 
use  of  an  explosive  source  essentially  requires  the  use  of  a  steel 
model,  and  steel  is  more  difficult  to  fabricate  into  models  than 
are  softer  materials  such  as  Plexiglas,  epoxy,  etc.  Furthermore, 
it  is  conceivable  that  the  force  of  the  explosive  charge  could 
affect  the  model  in  such  a  way  that  the  elastic  properties  of  the 
composite  layer  change  as  a  result  of  each  shot.  During  the  Rayleigh 
wave  work  special  precautions  were  taken  to  minimize  destruction  of 
the  model  by  the  source.  A  small  metallic  shield  was  placed  between 
the  source  and  the  model,  but  even  with  this  precaution ^ only  a 
limited  number  of  shots  could  be  fired  on  each  model. 
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The  second  reason  for  using  piezoelectric  devices  is 
to  have  a  repeatable  source.  With  a  pulsed  source  we  can  directly 
observe  the  seismic  records  on  the  oscilloscope  screen  and  can 
have  a  continuous  monitoring  of  the  wave  form  as  the  source  and 
detector  are  moved.  An  explosive  source  is  not  so  reproducible  as 
we  would  like  and  the  seismic  recoids  must  be  photographed  on  the 
oscilloscope  screen  to  be  observed.  In  addition,  the  piezoelectric 
sources  give  better  amplitude  reproducibility  than  the  explosive 
source.  The  principal  disadvantage  of  the  piezoelectric  sources 
and  receivers  is  that  they  do  not  give  the  sharp  wide -band  records 
of  the  explosive  source  and  that  the  resulting  seismograms  require 
heavier  filtering  than  the  explosive  records, 

Pigure  11  shows  a  block  diagram  of  the  seismic  modeling 
equipment  using  the  piezoelectric  equipment.  The  source  and 
receiver  are  identical  rectangular  pieces  of  lead  zirconium 
titanate  (PZT4)  manufactured  by  the  Cleavite  Corporation.  The 
source  is  activated  by  a  voltage  pulse  from  a  1,000  volt  pulser. 

The  wave  shape  from  the  pulser  is  variable  from  a  spike  of  1  micro¬ 
second  width  to  a  "fish  tail"  pulse  wirh  a  decay  time  of  5  milli¬ 
seconds.  The  voltage  pulse  most  used  in  our  work  was  a  fish  tail 
pulse  with  a  decay  time  of  one  millisecond.  The  basic  pulser 
unit  was  constructed  at  Calresearch,  and  the  circuit  diagram  is 
shown  in  Figure  12.  We  have,  in  addition,  two  commercially  made 
pulsing  units  with  an  output  of  30  volts.  These  units  are  useful 
for  records  made  at  small  offsets  from  the  source.  The  onset  of 
the  source  pulse  is  put  on  the  seismogram  by  a  Hewlett  Packer 


PULSER 
1000  VOLTS 


FIGURE  11 


SEISMIC  MODELING  PULSER 


-  58  - 


Digital  Delay  Generator,  Model  218a.  When  the  pulser  fires,  a 
low  voltage  signal  is  fed  into  the  digital  delay  generator.  At 
a  specified  number  of  microseconds  after  receiving  the  Input,  the 
generator  puts  out  a  signal  which  is  fed  directly  to  the  oscillo¬ 
scope  trace.  The  sweep  of  the  oscilloscope  is  triggered  by  the 
firing  of  the  pulser,  but  the  onset  of  the  sweep  is  not  reliable 
enough  to  be  used  as  a  time  reference.  The  signal  from  the 
receiver  is  first  put  through  a  buffer  amplifier  which  performs  the 
dual  function  of  matching  impedances  in  the  receiver  circuit  and 
providing  a  variable  gain  of  either  5  or  10.  The  circuit  diagram 
of  this  amplifier  is  shown  in  Figure  13.  After  amplification  the 
signal  is  fed  into  a  variable  band -pass  SKL  filter.  It  then  goes 
directly  into  a  Tektronix  5^^  oscilloscope.  The  digital  delay 
generator  is  adjusted  until  the  oscilloscope  trace  shows  a  time 
mark  some  5  microseconds  prior  to  the  onset  of  the  wave  form 
recorded  by  the  receiver. 

Figure  l4  is  a  photograph  of  the  modeling  equipment  when 
in  operation.  The  seismic  model  is  located  in  the  foreground.  It 
lies  flat  on  a  large  table,  wheheas  the  model  was  vertical  in  the 
explosive  source  work.  There  is  a  cushioning  layer  of  polyurethane 
foam  between  the  model  and  the  table  for  the  purpose  of  damping  out 
flexural  waves  in  the  model.  A  further  modifj. cation  to  the  model 
is  that  the  milling  of  the  steel  sheet  is  no  longer  used  in  the 
construction  of  the  composite  layers.  It  was  I'ound  that  if  a 
copper  strip  was  cemented  directly  to  the  steel  sheet,  a  low 
velocity,  composite  layer  wan  forir.ed  with  properties  essentially 
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identical  to  the  layer  in  which  a  section  of  the  model  was  first 
milled  out  before  the  addition  of  the  copper  strip.  The  model 
shown  in  the  photograph  consists  of  a  steel  sheet  8  feet  long, 

4  feet  wide,  and  .030  Inches  thick.  A  copper  strip  8  feet  long, 

1  inch  wide,  and  .030  inches  thick  was  attached  along  one  edge. 

This  technique  makes  it  relatively  easy  to  construct  either  layers 
or  more  irregularly  shaped  regions  of  different  elastic  properties 
for  the  seismic  modeling  investigations,  A  more  rig*d  epoxy  than 
was  used  in  the  explosive  work  was  found  to  be  best  in  the  con¬ 
struction  of  the  models  for  the  piezoelectric  work.  The  semi¬ 
rigid  epoxy  used  previously  caused  an  undesirable  attenuation  of 
energy  in  the  model.  The  epoxy  giving  the  best  results  was  a  rigid 
Hysol  epoxy,  base  A94899  and  hardener  H23044. 

In  Figure  14,  the  rack  immediately  behind  the  model  con¬ 
tains  the  filter,  the  digital  delay  generator,  and  two  different 
pulsing  units.  The  recording  oscilloscope,  with  camera  attached, 
is  adjacent  to  the  instrument  rack.  The  second  oscilloscope  is 
used  either  as  an  auxiliary  amplifier  or  to  monitor  the  unflltered 
wave  form.  The  source  is  located  In  the  holder  on  the  left  edge 
of  the  model,  while  the  receiver  is  on  the  long  rod  which  is 
between  the  instrument  rack  and  recording  oscilloscope  in  the 
photograph.  The  buffer  amplifier  is  located  just  beneath  the 
receiver. 

Figure  15  gives  two  examples  of  the  wave  forms  recorded 
with  the  equipment  shown  in»  the  previous  photograph.  The  first 
illustration  shows  the  direct  wave  and  the  first  reflected  arriv  .1 
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for  a  path  through  a  brass  sheet.  The  travel  path  is  30  centimeters 
long.  The  second  illustration  shows  a  refraction  record  for  a 
single  layer  at  an  offset  of  100  centimeters.  In  both  Illustrations 
there  is  evidence  that  resonance  in  the  source  and  receiver  con¬ 
tribute  to  the  wave  form  observed  on  the  seismic  record.  The 
source  and  receiver  have  a  primary  resonance  of  250  kc.  Filtering 
can  be  used  to  reduce  the  energy  at  this  frequency.  Frequency 
filtering  cannot,  however,  entirely  eliminate  resonance  phenomena, 
and  the  wave  shapes  are  never  so  sharp  as  in  the  work  done  with 
explosive  sources  and  strain  gage  receivers. 

The  primary  purpose  in  setting  up  refraction  modeling 
with  the  piezoelectric  sources  was  to  investigate  the  behavior 
of  the  refracted  arrival  through  a  high-velocity  layer  in  a  solid 
medium.  Theoretical  work  on  this  problem  is  summarized  in  another 
section  of  this  report''.  The  case  we  examined  by  modeling  tech¬ 
niques  was  that  of  a  high-velocity  layer  lying  between  two  lower 
velocity  layers,  with  these  three  laj/ers  overlying  a  high-velocity 
half -space.  We  examined  the  behavior  of  the  refracted  arrival 
through  the  high-velocity  layer  as  a  function  of  the  thicknesses 
and  elastic  properties  of  the  adjacent  layers.  We  accordingly 
began  with  a  study  of  the  propagation  through  a  system  consisting 
of  a  single  low-velocity  surface  layer  overlying  a  high-velocity 
half -space.  At  the  conclusion  of  this  work  two  additional  layers 
were  placed  Immediately  beneath  the  Icw-veloclty  layer,  and  the 
three-layer  profile  was  compared  with  t.hat  for  the  single  layer. 
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Figure  l6  shows  records  taken  every  5  centimeters  for 
receiver  offsets  of  from  5  to  200  centimeters  along  the  single 
layer  model.  The  direct  arrival  through  the  layer,  the  refracted 
arrival  through  the  half -space,  and  the  Rayleigh  wave  arrival  are 
all  clearly  defined.  Figure  17  gives  the  travel-time  curve  for  the 
direct  and  refracted  arrivals.  The  arrival  times  of  the  two  phases 
are  particalarD.y  well  defined.  There  is  very  little  scatter  in 
the  data.  The  velocities  deduced  from  the  travel-time  curve  are 
in  agreement  with  the  wave  velocity  measured  directly  through  both 
the  composite  surface  layer  and  the  steel  sheet.  The  layer  thickness 
inferred  from  the  travel  time  curves  is  23  mm.  This  is  approximately 
10  per  cent  lower  than  the  measured  tnickness  of  25.4  mm. 

In  addition  to  a  knowledge  of  the  arrival  times  of  the 
phases,  we  wish  to  determine  the  amplitude  behavior  of  these  arrivals. 
Data  on  amplitudes  are  not  nearly  so  reliable  as  data  on  arrival 
times.  The  principal  reason  for  this  uncertainty  is  the  effect  of 
the  variable  coup.llng  of  the  receiver  to  the  model.  Receiver 
positioning  is  all  important  in  modeling  work.  This  results  from 
the  relatively  large  size  of  the  receiver  when  compared  to  the 
other  dimensions  in  the  system.  Even  though  our  receivers  are  only 
two  millimeters  square,  they  correspond,  when  scaled  up  to  the 
real -earth  case,  to  seismometers  130  feet  square.  The  loading 
effect  of  such  seismometers  is  obviously  more  important  than  that 
of  a  short-period  Benioff  seismometer  in  actual  earthquake  recording. 
In  order  to  have  some  knowledge  of  the  relative  coupling  factor  at 
different  receiver  locations,  an  auxiliary  calibration  source  is 
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used  in  the  modeling  work.  This  source  is  placed  on  the  edge  of 
the  model  opposite  to  that  along  which  the  refraction  profile  is 
made  and  is  located  midway  between  the  source  and  the  receiver. 
Records  are  then  made  of  the  wave  form  of  the  calibration  source 
as  recorded  by  the  profile  receiver  and  by  the  profile  source 
used  as  a  receiver.  A  diagram  for  the  positioning  of  the  calibra¬ 
tion  source  is  shown  below. 


The  source  used  for  the  seismic  profile  is  fixed  at  S.  The  profile 
receiver  is  moved  along  the  various  positions,  Rj^ .  The  calibra¬ 
tion  source  is  moved  along  the  opposite  side  of  the  model  and  at 
each  receiver  position,  ,  is  located  at  a  position,  ,  which 
is  equidistant  from  S  and  R;^ .  The  object  is  to  calibrate  the 
receiver  at  each  position  along  the  profile. 

At  any  receiver  point,  R-^,  the  amplitude,  A^,  recorded 
on  the  oscilloscope  trace  is  given  by  ,  where  is  the 

"true"  wave  amplitude  immediately  below  the  receiver  position  and 
Mi  is  the  magnification  of  this  true  motion  by  a  combination  of 

V 
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the  receiver  plus  the  electronics.  Th.  assumption  is  made  that 
the  calibration  source  is  nondlrectional  so  that  for  any  Cj^ 
and  are  equal.  Wg.  is  the  wave  amplitude  at  S  due  to 
a  calibration  source  at  is  the  wave  amplitude  at  R^due 

to  the  source  at  .  V/e  wish  to  determine  the  ratio  Mp  . 

We  can  measure  the  initial  amplitude  of  the  two  calibration 
records  recorded  at  S  and  and  can  form  the  ratio  of  these 
amplitudes 

y  = 

This  gives  then  that  M^/M^  =  assumptions  that 

Me  =  Mo  #  since  the  source  is  fixed  at  S,  and  that  Wo  =  Wo 

o, 

because  the  source  is  nondlrectional.  Using  the  records  from  the 
calibration  source  to  normalize  the  refraction  profile  ampli¬ 
tudes,  it  is  possible  to  get  a  reasonably  reliable  measure  of  the 
amplitudes  at  the  various  receiver  positions. 

Figure  l8  shows  the  normalized  refracted  wave  amplitudes 
for  the  single  layer  model.  These  amplitudes  do  not  fit  the 
theoretical  predictions  (Heelan,  195^)  of  an  attenuation  as  the 
inverse  3/2  power  of  receiver  offset,  but  show  a  more  rapid  attenua 
tion  with  distance.  An  attenuation  varying  as  the  -1.75  power  of 
receiver  offset  fits  the  observed  data.  Three  profiles  were  run 
along  the  single  layer  model  and  all  gave  similar  attenuation 
coefficients  for  the  refracted  wave. 


AMPLITUDE 


FIGURE  18 

REFRACTED  WAVE  AMPLITUDE-ONE  LAYER  MODEL 


LE  44-542 
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Figure  19  shows  the  profile  along  a  model  consisting  of 
three  layers  overlying  a  half -space.  The  direct  and  refracted 
wave  velocities  are  the  same  as  for  the  single  layer  case.  At  the 
time  this  profile  was  taken,  the  receiver  positions  unfortunately 
were  not  calihrated,  and  the  amplitude  data  are  not  sufficiently 
accurate  for  comparison  with  the  single  layer  case. 


RECEIVER  OFFSET 
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IV.  LONG-TIME  RESPONSE  PREDICTED  BY 
EXACT  ELASTIC  RAY  TliEORY 

A.  Introduction 

Cagniard  (1939)  analyzed  In  considerable  detail  the 
response  produced  by  the  reflection  of  a  spherically  symmetric 
compressional  wave  at  a  plane  Interface  between  two  solid  half¬ 
spaces.  The  total  response  is  the  sum  of  a  compressional  component 
and  a  shear  component,.  For  a  step  function  input  in  particle 
velocity,  both  components  diverge  in  the  long-time  limit.  Each 
component  depends  on  an  integral  which  goes  to  zero  in  the  long¬ 
time  limit  and  on  a  polynomial  in  time  which  contains  a  linear 
term  and  a  cubic  term.  When  the  compressional  and  shear  components 
are  added  together,  the  cubic  terms  cancel  one  another  but  the 
linear  terms  remain.  The  linear  divergence  results  from  using  an 
input  function  (a  step)  which  is  not  permissible  on  physical 
grounds.  When  a  physically  realizable  function  is  used,  the  total 
response  goes  to  zero  in  the  long-time  limit  even  though  the  com¬ 
ponent  parts  of  the  total  response  diverge.  Cagniard ‘s  analysis 
is  here  extended  to  rays  which  are  multiply  reflected  in  a  multi¬ 
interface  system,  and  a  prescription  is  presented  for  arranging 
the  rays  In  groups.  Each  group  response  function  goes  to  zero  in 
the  long-time  limit  even  though  its  component  parts  diverge.  The 
prescription  facilitates  the  computatxcn  of  the  exact  response 
function  over  an  extended  time  Interval  in  which  many  rays  contri¬ 


bute  . 
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The  vertical  component  of  particle  velocity  in  a  spheri¬ 


cally  symmetric  compressional  wave  is 


) 


T-t- 


(1) 


The  plane  Z  =  0  is  at  a  depth  H  below  the  source  and  at  a  depth 
Z  below  the  receiver  (Figure  1).  is  the  source-receiver 

separation.  V  is  the  compressional  wave  velocity  and  b  is  the 
density.  Subscript  one  refers  to  medium  one.  gl  is  the  radius 
of  the  spherical  cavity  and  is  the  peak  pressure  applied  to  the 
cavity  wall-  The  choice  of  the  source  function,  ^  ,  is  not  com¬ 
pletely  arbitrary.,  It  must  be  chosen  so  that  each  point  in  the 
medium  returns  to  its  original  state  a  finite  time  after  the 
arrival  of  the  disturbance  at  that  point.  This  implies  that  the 
static  component  of  such  quantities  as  the  stress  and  strain  must 
be  finite.  The  Laplace  transform  of  the  radial  stress  produced 
by  a  spherically  symmetric  compressional  wave  is 


W  r,  HHSPoj  ^ 


4'  is  the  Laplace  transform  of  4"  »  S  is  Laplace  transform 

variable  and  'Zr  is  the  shear  velocity.  In  order  for  the  static 

component  of  stress  to  be  finite 

limit  )<  .=  2.  - 

S-^0  ‘  ^ 

As  a  consequence  of  (2),  )f(T)  niust  hive  at  least  two  axis 
crossings  and 


(2) 


oO 


J 


[  (C.4.C2>U/;)f(t-^)4^ 


(3) 


—  od 


for  arbitrary  and 


^  j{ 

I  H  " 


An  axially  syiTimetric  radiation  field  in  a  layered  system 
can  be  expanded  in  an  infinite  series.  Each  term  can  be  derived 
directly  from  the  integral  representation  for  the  Laplace  transform 
of  the  source  radiation  field  by  using  the  method  of  generalized 
reflection  and  transmission  coefficients  (Spencer,  i960).  This 
method  provides  a  simple  prescription  for  transforming  the  integral 
representation  for  the  source  radiation  field  into  an  integral 
representation  for  the  response  function  associated  with  any 
generalized  transmission  path.  The  quantities  which  are  required 
to  specify  a  generalized  transmission  path  contain  all  the  informa¬ 
tion  required  by  the  prescription  to  construct  the  associated 
Integral  representation.  A  particular’ generalized  path  is  com¬ 
pletely  determined  by  specifying;  (a)  the  total  vertical  distance 
)  traveled  in  each  layer  in  each  mode  ( compresslonal  and  shear) 
and  (b)  the  sequence  in  which  the  layers  are  traversed. 

The  integral  representation  in  cylindrical  coordinates 


of  the  Laplace  transform  of  (1)  is 


pis  the  radial  distance  from  the  axis  of  symmetry  (Z-axis). 

is  the  zero  order  Besses  "^unction  of  the  first  kiiid.  The  pre¬ 


scription  gives  the  Laplace  transform  of  the  response  associated 


with  each  generalized  path  in  the  form 


w 
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M  is  twice  the  number  of  layers  traversed  by  the  disturbance.  The 
are  the  compresslonal  and  shear  wave  velocities  in  the  differ¬ 
ent  layers.  A( ^  ,  S  )  represents  the  product  of  the  generalized 
reflection  and  transmission  coefficients. 


B .  Long-Time  Behavior 

Cagniard  has  developed  an  elegant  method  for  inverting 
Laplace  transforms  of  the  form  given  in  (4).  S  is  treated  as  a 
positive  real  variable.  This  restriction  on  S  is  a  distinctive 
feature  of  Cagniard *s  theory.  The  transformation  A=SfyV|^ 
in  (4)  gives 

A(f  ,  K(-f  >S)=5|  )  •  (5) 

1-1 


Because  S  is  real  this  transfonr.atlon  is  nothing  more  than  a 
change  in  scale,  yet  it  serves  to  vastly  sliipllfy  the  S-dependence 
of  the  Intefjrand  and  to  express  that  dependence  in  explicit  form. 

^  is  dimensionless.  and  refer  to  a  characteristic 

dimension  and  velocity  respectively. 

The  branch  points  of )  and  ^  )  are  associated 
with  the  square  roots  (  The  branch  points  of-pCf  ) 

are  labeled  t  and  those  of  l(f )  are  labeled 

i  (l«c  1,-^M)  .  The  condlt-’on  •  '  •  '  is  always 

satisfied.  The  integrand  is  made  single- valued  by  cutting  the 
^  -plane  along  the  imaginary  axis  between  each  branch  point  and 
its  c'^mplex  conjugate  (Figure  2).  We  confine  our  attention  to 
the  sneet  of  the  Riemann  su’^face  on  which  all  the  square  roots 
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are  positive  real  on  the  positive  real  £  -axis .  Cagniard  demon¬ 
strated  that  the  branch  points  are  the  only  singularities  except 
in  the  very  special  case  where  an  interface  can  support  a  Stoneley 
wave.  For  each  interface  which  can  support  a  Stoneley  wave,  (f*  ) 
has  two  pure  imaginary  conjugate  poles  CQ,^  ).  The  order  of 

the  poles  associated  with  a  particular  Interface  is  determined 
by  the  number  of  times  the  disturbance  Interacts  with  (is 
reflected  from  or  transmitted  across)  that  Interface. 


Applying  Cagniard ’s  method  to  (4)  gives 


F{T)  is  the  vertical  component  of  the  particle  velocity  produced 
by  a  step  variation  in  the  source  function.  When  (the 

travel  time  via  the  least -time  reflection  path),  F(T )  can  be 


expressed  in  the  form 


where 


FJT)=At?f  , 


I  indicates  that  only  the  real  part  of  the  quantity  in  paren¬ 
thesis  is  required.  The  contour,  B,  wraps  around  the  part  of  the 


branch  cut  between  -i4j  and  -  in  a  clockwise  sense  (Figure  2.) 
Any  poles  below  -  ^  are  enclosed  in  a  clockwise  sense  and  are 

included  in  B.  In  the  long-time  limit  F^  decays  like  H  *  . 
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The  I  contour  is  a  seni-circle  in  the  lower  half  plane  with 
center  at  the  origin.  The  contour  is  traversed  in  a  counterclock¬ 
wise  sense.  The  value  approached  by  the  Integral  in  (8)  in  the 
limit  as  the  radius  of  the  semi-circle  goes  to  infinity  determines 


Fj  (T). 

A  careful  examination  of 

lotut  I  _ L  y  f  ^nCt) 

n-0  \  ^ 


wnere 


reveals  that 


(T  )  and  p2n+|  ^  polynomials  in  even  and  odd  powers  of  > 

respectively. 

The  form  of  -P  (  §*  )  depends  on  whether  the  generalized 
transmission  path  is  degenerate.  Degeneracy  occurs  when  one  or 
more  of  the  layers  are  traversed  in  both  the  compressional  and 
shear  modes.  Then  there  are  two  or  more  generalized  paths  for 
which  the  (and  consequently  the  functions  ))  are  identical 
but  the-P{.§  )  are  in  general  different.  We  can  group  all  degenerate 
paths  together  by  modifying  the  form  of  4*  (i*  ) •  The  response 
functions  associated  with  degenerate  paths  all  have  identical  onset 
times.  When  the  path  is  not  degenerate^  4' (i*  )  is  expressed  by  a 


,  7ft  _ 

f 

product  of  the  generalized  coefficients.  When  the  path  is  degen¬ 
erate,  -P(f  )  is  expressed  as  a  sum  of  products  of  the  generalized 
coefficients . 

The  behavior  of  -p  ( f  )  foT’  large  f  can  be  derived  from 
the  expressions  for  the  generalized  reflection  and  transmission 
coefficients.  Each  of  the  four  generalized  reflection  coefficients 
has  the  form 

lOTit  R(£)=r£"-t-rfVrrV  •  ■  •  • 

If  the  incident  wave  approaches  the  interface  through  a  solid 
(2Jr¥’0  ),  cannot  vanish  in  any  of  the  reflection  coefficients. 

If  the  incident  wave  approaches  the  Interface  through  a  fluid,  there 
is  only  one  reflection  coefficient  (  i^pp.  )*  and  .  If  both 

media  are  fluids  and  there  is  no  density  contrast,  r^=rO  •  The 
four  generalized  transmission  coefficients  have  the  form 

Limit ■ 

may  vanish  under  certain  conditions. 

If  the  trar.- ..xsslon  path  is  not  degenerate,  -F(S  )  is 
Just  the  product  of  the  generalized  coefficients.  In  what  follows 
we  assume  that  there  are  no  fluid  layers.  Let  nr)  be  the  number  of 
reflection  coefficients  in  the  product.  If  a  reflection  coefficient 
is  raised  to  the  Jfi,  power,  it  is  counted  Jk  times  in  determining  rf)  • 


♦The  letter  subscripts  indicate  the  mode  of  propagation  before  and 
after  transmission  (P  for  compressionai  and  S  for  shear).  The  number 
subscripts  Indicate  the  layers  in  which  the  disturbance  propagates 
(the  layer  is  above  the  interface). 
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Let  ^  be  the  number  of  transmission  coefficients  for  which  the 

term  vanishes.  if>U'jr«.o  ,  in  a  transmission  coefficient  which 
is  raised  to  the  X  ^  power,  that  coefficient  is  counted  X  times  in 
determining-^.  If  the  transmission  path  is  degenerate, -f  )  will 
oe  given  by  a  sum  in  which  each  term  consists  of  a  product  of  the 
generalized  coefficients.  Both  the  number  of  reflection  coeffi¬ 
cients  and  the  number  of  transmission  coefficients  will  be  the  same 
for  each  term  in  the  sum.  Consequently,  the  same  ■'’alue  of  fT\  must 
be  associated  with  each  term  in  the  sum.  However,  even  though  the 
number  of  transmission  coefficients  is  the  same  in  each  term,  the 
particular  coefficients  need  not  all  be  the  same.  Consequently, 
the  value  of-^  need  not  be  the  same  for  every  term.  Let  be  the 
smallest  of  the  ^  values.  Then 

regardless  of  whether  the  path  is  degenerate  or  not.  The  are 
functions  of  the  density  and  velocity  ratios  across  the  interfaces 
intersected  by  the  generalized  transmission  path. 

For  large  £  the  integrand  in  (8)  can  be  expressed  as 
an  infinite  series  in  powers  of  ^  .  When  integrated  around  the 
semi -circle,  even  powers  of  ^  give  pure  imaginary  values  which 
contribute  nothing  to  F^(T)-  Except  for  the  term,  odd 

powers  of  ^  vanish  when  Integrated  around  the  semi -circle. 
Integration  of  the  term  gives 

F^(T)-0  J  ^ 

m^o+1 

When  inricfjj,,  vanishes  because  there  is  no  term. 
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Each  of  the  P  «  ( T )  contain  only  odi  powers  of  T  . 

There  are  A  terras  In  ^  (-T  )  and  T  is  the  highest  order  terra. 
Therefore,  F  (T  )  can  also  be  expressed  in  the  fcrra 


F^(T)=  2^  .  (10) 

1-0 

The  determination  of  the  Is  straightforward  but  tedious  even 
for  simple  transmission  paths. 

The  long-time  response  is  determined  by  a  polynomial  of 
degree  in  odd  powers  of  T  .  Suppose  ■^o"0  .  For  direct 

transmission  paths,  lflf)=0  and  F^  contains  only  the  linear  term. 

For  primaries,  rfWl  and  F^  contains  both  a  linear  and  a  cubic  term. 
For  an  order  multiple,  F^  contains  all  the  odd  powers  of  T 
through  T  .  This  divergent  behavior  is  a  characteristic  feature 
of  the  response  functions  for  the  individual  transmission  paths. 


C.  Convergent  Groups 

The  existence  of  the  divergent  tails  can  be  extremely 
troublesome  in  the  numerical  evaluation  of  the  exact  response  over 
an  extended  time  interval  in  which  higher  order  multiples  contribute. 
For  example,  in  using  a  digital  computer  to  evaluate  and  sura  the 
individual  response  functions,  an  Intermediate  stage  of  calculation 
may  be  reached  in  which  the  magnitude  of  the  partial  sum  exceeds 
the  magnitude  of  the  correct  total  sura  by  a  number  which  exceeds 
the  number  of  significant  figures  carried  by  the  computer.  When 
this  occurs  the  final  computed  sum  is  nonsense.  Pekerls  and 
Longman  (1958)  did  not  encounter  this  difficulty  for  the  case  of 
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a  fluid  layer  over  a  fluid  half -space  because  =  0  In  the 
expansion  for  the  reflection  coefficient.  Consequently,  the 
constant  term  Is  the  leading  term  In-P  (£  )  for  all  paths.  It 

follows  that  the  response  produced  by  an  Input  step  dl\>  -‘ges 

linearly  and  the  response  produced  by  an  Impulse  approaches  a  con¬ 
stant  for  each  path.  For  the  all  solid,  multi-layer  systems,  we 

can  circumvent  this  problem  by  collecting  the  ray  paths  Into  groups 
Equation  (6)  shows  that  the  long-time  response  produced  by  a 
particular  source  function  Is  obtained  by  convolvjng  with  d^/dt 
Equation  (3)  shows  that  for  a  physically  realizable  source  function 
dFj.  /dt  must  not  diverge  more  rapidly  than  i:  ,  otherwise  the 
medium  Is  left  with  a  residual  particle  velocity.  This  suggests 
that  the  rays  can  be  grouped  together  In  such  a  way  that  the  cubic 
and  higher  order  terms  vanish  In  the  group  response  to  a  step 
function  Input, 

By  analyzing  several  cases  we  have  found  that  the 
generalized  transmission  paths  v/hlch  belong  to  a  particular  group 
have  the  following  characteristics  In  common: 

(a)  the  Interfaces  at  which  the  reflections  take  place, 
(B)  the  sv=iquence  In  which  the  reflecting  Interfaces 
are  encountered. 


(C)  the  generalized  path  to  the  first  reflector. 


Let  P  be  the  step  function  response  associated  with  the 

U 

generalized  transmission  path.  The  response  function  obtained  by 


summing  all  the  which  have  (A),  (B),  and  (C)  In  common  varies 
linearly  In  the  long-time  limit. 
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Figure  3  shows  four  primaries  which  are  reflected  from 
the  (l+i  interface:  P.  ^  S-  P.  P.  .  ,  P.  S  P  S.  ,  P.  S.  a  P  , 
and  P^i  a  j  .  These  primaries  have  a  common  transmission  path 
to  the  reflector  and  form  a  group.  When  we  say  that  the  paths  are 
common  down  to  the  reflector,  we  mean  that  the  mode  of  transmission 
across  each  layer  between  the  source  and  the  reflector  is  the  same. 
This  does  net  mean  that  the  least-time  reflection  paths  are  the 
same  -  in  fact,  they  are  different.  This  group  does  not  include 
all  the  primaries  from  the  (t-tl  )  interface.  Actually  there  are 
sixteen  primaries  which  can  be  arranged  In  four  groups  of  four 
each.  All  members  of  each  group  have  common  paths  down  to  the 
reflector.  Each  group  response  function  varies  linearly  in  the 
long-time  limit.  Consequently,  he  total  primary  response  func' 
tlon  varies  linearly  in  the  long-time  limit. 

Now  consider  the  general  case  wher  ■  the  paths  suffer  fT) 

reflections.  Let  the  reflectors  be  designated  by  . .  . .  ,ffl  ) 

and  let  the  order  in  which  the  R  's  are  arranged  indicate  the 

i\ 

sequence  in  which  the  reflectors  are  encountered.  If  the  dis¬ 
turbance  is  reflected  from  the  same  Interface  more  than  once,  the 
designation  for  that  Interface  will  appear  more  than  once  in  the 

sequence.  The  total  number  of  pat'^s  which  go  through  a  parti- 
i\ 

cular  reflection  sequence  and  have  a  common  path  to  the  first 
reflector  is  2  •  ,  Nj  is  the  number  of  interfaces  which  all 

paths  encounter  in  going  from  the  Rj  reflector  to  the  receiver. 

The  highest  order  term,  in  the  Individual  is  T  .  The  high¬ 
est  order  term  rU  oained  by  summing  all  the  which  have  a  common 
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path  tc  the  last  reflector  Is  T  '  .  The  highest  order  term 

obtained  by  summing  all  the  which  have  a  common  path  to  the 
next  to  last  reflector  is  .  Each  time  we  reduce  by  one 

the  number  of  reflectors  on  the  common  part  of  the  path,  we  reduce 
the  exponent  of  the  highest  order  term  in  the  resultant  polynomial 
by  two  and  increase  the  number  of  contributing  paths  from  < 

li"! 

to  2  .  The  highest  order  term  in  the  function  obtained  by 

summing  all  the  whicn  have  a  common  path  to  is  O'  . 

In  computing  the  response  of  a  layered  system  we  first 
determine  those  generalized  paths  for  which  0^  lies  within  the 
time  Interval  of  interest.  These  paths  are  arranged  in  groups 
according  to  (A),  (B),  and  (C).  In  each  group  we  determine  the 
0^  which  is  largest  (0^ )  and  evaluate  the  individual  F^  in  the 
interval  ,  where  T  is  the  duration  of  the  source 

function.  There  is  no  need  to  calculate  the  F^  for  larger  times 
because  their  sum  varies  linearly  and  vanishes  when  convolved  with 
a  physically  realizable  source  function.  The  conditions  (A), 

(B),  and  (C)  determ.ine  the  minimum  number  of  rays  which  when  added 
together  have  the  desired  long-time  behavior.  By  grouping  the 
individual  rays  according  tc  (A),  (B),  and  (C),  we  can  construct 
convergent  functions  from  divergent  ones.  The  total  response  pro¬ 
duced  by  a  physically  realizable  source  function  can  be  expressed 
in  terms  of  functions  which  decay  in  the  long-time  limit  even 
though  the  component  parts  of  these  functions  diverge. 


/Ok  • 


^^11  - 

I 


^  -  PLANE 


FIGURE  2 


LE44-828 


Source  - 


r; 

‘^^-1 

r, 

^S^-l 

'  py 

*- 1 

i-\  Medium 

r-  ' 

i 

i  1  Medium 

FIGURE  3 


Receiver 


LE  42-618 


-  87  - 


V.  HIGH-FREQUENCY  ELASTIC  WAVE  THEORY 


A .  Introduction 

Because  of  the  mathematical  complexity  of  even  the  simplest 
problems  In  elastic  wave  theory,  considerable  work  has  been  directed 
toward  obtaining  approximations  to  the  exact  theory.  Two  of  the 
most  useful  approximations  are  the  normal  mode  theory  (valid  at 
large  ranges)  and  geometric  ray  theory  (valid  at  high  frequencies). 

In  using  the  approximations  we  generally  do  not  know  what  the 
error  Is  -  'll!  we  know  is  that  In  some  limit  the  error  is  negligible 
and  that  the  limit  itself  depends  on  how  we  choose  to  measure  the 
error.  The  limit  Is  usually  estimated  by  determining  the  point 
at  which  some  prediction  of  the  approximate  theory  begins  to  depart 
from  what  seems  physically  reasonable.  Our  objective  is  to  develop 
the  high  frequency  theory  for  a  layered  medium  and  to  point  out 
some  of  its  shortcomings. 

B .  The  Source  Function 

We  consider  the  response  produced  by  a  source  which 
radiates  a  spherically  symmetric  compressional  wave.  In  an  infinite 
homogeneous  medium,  the  vertical  component  of  the  particle  velocity 
is 


where 

T--  (f\cz-wfy^  * 
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p  is  the  source-receiver  separation.  The  plane  Z—0  is  at  a  depth 
H  below  the  source  and  Z  below  the  receiver.  is  the  peak  pres¬ 
sure  applied  to  the  cavity  wall  and  a  is  the  radius  of  the  cavity. 

and  refer,  respectively,  to  the  compress ional  velocity, 
shear  velocity  and  density  in  medium  one.  The  source  function, 

>f,  merely  expands  or  contracts  with  "t^  but  does  not  change  shape. 
The  amplitude  spectrum  of  is  assumed  to  peak  at  a  dominant 
frequency  .  At  distances  which  are  large  compared  with 

the  dominant  wavelength,  the  second  term  in  (1)  becomes  negligible 
with  respect  to  the  first  term,  the  disturbance  propagates  without 
change  of  shape,  and  ^ (T )  determines  the  time  variation.  In 
effect,  )f ( O'  )  is  a  high-frequency  apnroxlmatlon  to  the  response 
of  an  infinite  medium. 

To  determine  the  range  of  validity  of  the  high-frequency 
approximation  we  first  demonstrate  that  and  the  Integral  of  )( 
take  on  non-zero  values  over  the  same  Interval.  X  must  be  chosen 
so  that  each  point  returns  to  its  original  state  a  finite  time 
after  the  arrival  of  the  disturbance  at  that  point.  This  means 
that  the  static  component  of  the  stress  and  strain  must  be  finite. 
The  Laplace  transform  of  the  radial  stress  in  a  spherically 
symmetric  compressional  wave  is 


In  order  for  the  static  component  to  be  finite, 
Urnit  )(  -  AS^  ^^-2  ■ 


(2) 
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As  a  conseqjence  of  (2),  the  static  component  of  }(  must  vanish. 

This  insures  that  the  integral  of  )f  is  zero  outside  the  interval 
where  X  takes  on  non-zero  values.  Therefore,  the  validity  of  the 
high-frequency  approximation  depends  only  on  the  relative  amplitudes 
of  the  first  and  second  terms  in  (1).  If  the  two  time  functions 
had  been  resolved,  their  significance  would  depend  not  on  their 
relative  amplitudes  but  on  the  amplitude  of  the  background  (signal 
plus  noise)  which  is  coincident  in  time.  If  there  were  no  back¬ 
ground,  neither  term  could  be  neglected.  For  purposes  of  illus¬ 
tration  let 

TfCT)  —(Sin^^T  )  ( COS  2n)  ^  K  ( 1^2)  . 

The  first  factor  determines  the  shape  of  the  envelope.  K  deter¬ 
mines  the  number  of  oscillations  and  the  bandwidth.  The  condition 
Ki^  2  forces  the  low  frequency  behavior  to  satisfy  (2).  For  K  =  4, 
the  ratio  of  the  maximum  amplitudes  is  .l6“l:^V|yr  .  Therefore, 
for  a  dominant  wavelength  to  range  ratio  of  one  half,  the  second 
term  is  about  eight  per  cent  of  the  first  term.  This  provides  an 
estimate  of  the  upper  limit  on  the  dominant  wavelength  to  dimension 
ratio  for  which  the  high-frequency  approximation  is  valid. 

C.  Generalized  Ray  Theory 

An  axially  symmetric  radiation  field  in  a  layered  system 
can  be  expanded  in  a  series.  Each  term  can  be  derived  directly 
from  the  Integral  representation  for  the  source  radiation  field 
by  using  the  method  of  generalized  reflection  and  transmission 
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coefficients  (Spencer,  i960).  This  method  provides  a  prescription 
for  transforming  the  integral  representation  for  the  source  radia¬ 
tion  field  into  an  integral  representation  for  the  response  func¬ 
tion  associated  with  any  generalized  transmission  path.  The 
quantities  which  are  required  to  specify  a  particular  generalized 
transmission  path  contain  all  the  information  required  by  the  pre¬ 
scription  to  construct  the  associated  integral  representation. 

A  particular  generalized  path  is  completely  determined  by  specify¬ 
ing:  (a)  the  total  vertical  distance  traveled  in  each  layer  in  each 
mode  ( compress lonal  and  shear)  and  (b)  the  sequence  in  which  the 
layers  are  traversed. 

The  integral  representation  in  cylindrical  coordinates 


of  the  Laplace  transform  of  (l)  is 


where 


The  prescription  gives  the  Laplace  transform  of  the  response  associa¬ 


ted  with  each  generalized  path  in  the  form 

'  o 

where 


(3) 
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The  summation  extends  over  all  layers  traversed  by  the  disturbance. 
The  and  T|^  are  the  total  vertical  distances  traveled  In  the 
layer  In  the  compresslonal  and  shear  modes  respectively.  ) 

represents  the  product  of  the  generalized  reflection  and  trans¬ 
mission  coefficients. 

The  functions  A{\S  )  and  K(^jS  )  contain  square  roots 
of  the  form  (  refers  to  a  body  wave  velocity. 

To  make  the  integrand  single-valued,  we  cut  the  A-plane  along  a 
straight  line  between  the  branch  points.  S  is  considered  complex 
of  the  form  S  — (a);>0).  The  angle  which  the  cut  makes  with 
the  real  axis  and  the  argument  of  the  square  root  at  points  on  the 
cut  is  indicated  in  Figure  1. 

As  the  cut  approaches  the  path  of  integration 

(along  the  positive  real  A -axis)  from  below.  Let  S—LcO  ,  then 
the  path  of  integration  lies  above  the  cut.  The  change  of  scale 
A— transforms  (3)  into 


In  (5)  we  have  used  one  expression  for  both  the  compressional  and 
shear  terms.  may  refer  to  either  or  and  M  is  twice  the 
number  of  layers  traversed  by  the  disturbance.  Hp  and  V„  may  be 
associated  with  a  characteristic  dimension  and  velocity  respec¬ 
tively. 

As  ^7^/2  i  the  cut  approaches  the  path  of  integration 

from  above.  Let  S  =  -LCO  ,  then  the  path  of  integration  lies  below 
the  cut.  The  change  of  scale  transforms  (3)  into 


X(-^) 


r 


y 


BEIDW 


(6) 


The  only  difference  between  the  integrals  in  (4)  and  (6)  is  in 
the  position  of  the  cut  with  respect  to  the  path  of  integration. 

U  is  the  transform  of  a  real  time  function  therefore, 

(u(iooj)^—  U  (-^^)  )  . 

The  asterisk  designates  the  complex  conjugate,  yf  is  by  defi¬ 
nition  a  real  time  function,  hence 

X(-cw)  . 


T.he  relationship 


follows  directly  from  (5). 
real  implies 


Therefore,  the  condition  that  L)  be 


-f(0  )  .  0:^§^oC, 
ABOVE/  ^  BELOVV/  ^ 


(7) 
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This  result  Is  one  of  a  group  of  symmetry  relations  which  relate 
the  values  of-f(^  )  above  and  below  the  cut  and  at  points 
reflected  In  the  Imaginary  axis.  The  symmetry  relations  are 
valid  for  real  ^  and  follow  directly  from  the  expressions  for 
the  generalized  coefficients. 

The  spectrum  for  an  Impulsive  source  function  1  )  Is 

(■^XvR)n=i'/Y(Y)) 

wnere 

Yh)- [  ,  Y>o. 

ABO^/E 

A  knowledge  of  the  nature  and  position  of  the  singularities  of 

)  and^(f  )  Is  necessary  for  the  development  of  the  asymptotic 
expansion  and  to  a  proper  undei-s tending  of  the  results.  We  con¬ 

fine  our  attention  to  the  sheet  of  the  Rlemann  sur^’ace  on  which 
the  square  roots  (  are  positive  real  on  the  real 

axis  to  the  right  of  the  branch  points.  All  the  singularities  of 
•f(f  )  and  §  )  lie  on  the  real  axis.  The  branch  points  of  -f(^  ) 
are  designated  by  t  =  1,...,N)  and  those  of/^(£  )  by 

If  the  body  wave  velocities  in  the  layers  are  all 
different,  the  branch  points  of4-(^  )  form  a  set  which  includes 
all  the  branch  points  ofy|^(i'  )  plus  additional  ones.  )  con¬ 

tains  the  four  branch  points  associated  with  each  interface 
Intersected  by  the  generalized  transmission  path.  contains 

only  branch  points  which  are  associated  with  the  particular  mode 
of  propagation  In  the  layers.  Even  when  the  generalized  trans¬ 
mission  path  traverses  every  layer  In  both  modes, -r(^)  contains 


four  additional  branch  points  associated  with  the  layers  immedia¬ 
tely  above  the  -ppermost  reflector  and  immediately  below  the  lower¬ 
most  reflector.  Note  that  if  there  is  no  reflection  at  the  inter¬ 
faces  between  which  the  transmission  path  is  contained,  the 
generalized  transmission  path  cannot  traverse  the  outermost 
layers  in  both  modes.  Then,  even  if  all  the  other  layers  were 
traversed  in  boch  modes, )  would  still  contain  two  branch 
points  not  contained  ).  If  the  body  wave  velocities  are 

not  distinct,  special  cases  can  arise  where  the  branch  points  of 
-fh  )  andi(f  )  are  identical.  The  elation 

is  satisfied  in  ail  cases. 

The  poles  of  the  integrand,  when  they  exist,  are  associa¬ 
ted  with  zeros  of  the  demoninators  in  the  expressions  f)r  the 
generalized  coefficients.  Cagniard  (1939)  demonstrated  that  at 
an  interface  between  two  elastic  half-spaces  the  denominator  has 
no  zeros  on  the  permissible  sheet  of  the  Riemann  surface  except 
in  the  very  special  ca..e  where  the  interface  can  supi  .’u  a  Stoneley 
wa\e.  Then  the  denominator  has  the  two  real  zeros  = 

JJj[  Is  the  Stoneley  velocity  on  the  Interface.  Spencer  (1956) 
showed  that  a Stoneley  wave  is  always  excltea  at  a  liqu  'd-solid 
interface.  A  Rayleigh  wave  is  always  excited  at  a  solid-vacuum 
interface.  In  all  three  cases  the  zeros  are  real  and  are  situated 
symmetrically  Vvith  respect  to  the  origin.  is  always  less  than 
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the  smallest  body  wave  velocity  in  the  adjacent  media.  Therefore, 
the  pole,  ,  cannot  lie  on  any  of  the  four  cuts  associated 


with  the 


interface  and 


minimum  . 

In  deriving  the  asymptotic  expansion  of  Y(n^),  we  find 


that  only  the  case  ^  =  0  can  yield  terms  associated  with 

l~i 

the  poles.  This  situation  arises  when  the  souri.3  and  receiver 
are  situated  in  the  same  horizontal  plane  and  the  generalized 
transmission  path  is  a  direct  path.  If  the  source  and  receiver 
are  in  the  same  plane  but  not  on  an  interface,  the  direct  paths 
do  not  interact  with  any  interface.  The  response  function  is 
then  Just  the  time  variation  in  the  direct  wave,  a  function  which 
is  prjscribed  in  the  initial  statement  of  the  problem.  If  the 
source  and  receiver  are  on  the  same  interface,  the  direct  paths 
interact  with  that  interface  only  and  the  problem  reduces  to  the 
one  interface  case.  Y(^  )  becomes 


Y(^)  =  f f J.(rf I ) -f (f)dS  ,  = 0 . 


Clearly  the  asymptotic  expansions  for  high  frequency  and  large 
offset  are  identical. 

By  expressing  the  Bessel  function  in  terms  of  the 
Hankel  functions  of  the  first  and  second  kinds  and  using  Cauchy’s 
Integral  theorem,  we  can  demonstrate  that 


M 


ABOVE 
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The  asymptotic  expansion  of  the  branch  line  integral  can  be 
obtained  by  using  a  variation  of  the  method  of  stationary  pha^e 
discussed  by  Erdelyi  (1956).  The  high  frequency  behavior  is  deter 
mined  by  the  four  branch  points  (N  =  4  for  the  one  Interface 
problem).  This  is  to  be  expected  because  the  path  along  the  inter 
face  is  a  least-time  path  for  four  modes  of  propagation  (i.e.,  the 
compresslonal  and  shear  body  waves  in  the  adjacent  media).  The 
second  term  in  (10)  is  present  only  if  the  interface  can  support 
a  Rayleigh  or  Stoneley  wave .  The  purpose  in  discussing  a  problem 
which  has  been  extensively  treated  in  ■:he  literature  is  to  point 
out  the  very  special  situation  under  which  the  pole  and  every 
branch  point  of-p(^  )  contribute  to  the  high  frequency  asymptotic 
expansion. 

When  ^  0,  the  expession  which  corresponds  to 


(10),  is 


sine 

-sunO 


4^  1  rjs 


SUlG 


(C+) 

Q  is  the  angle  between  the  vertical  and  the  segment  of  the  least¬ 
time  reflection  or  transmission  path  which  is  traversed  with 

velocity  V  .  sinG  satisfies 
K 


T.  ^  sin0  y  _ 


0^  . 


(12) 
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The  mathematical  significance  of  sln0  can  be  seen  by  replacing 
the  Hankel  functions  by  their  asymptotic  representations.  In 
the  resultant  integrands,  appears  only  in  the  exponentials 
S  and  0.  ^ 


where 


£_  —-sin©  is  the  only  saddle  point  on  the  surface  | 

in  the  upper  half  plane.  ^  =  sin©  is  the  only  saddle  point  on 
the  surface  in  the  upper  half  plane.  The  contours 

C_  and  can  be  continuously  deformed  into  lines  slightly  above 
the  part  of  the  real  axis  to  the  right  of  the  saddle  points  with¬ 
out  enclosing  any  singularities  of  the  integrand. 

When  sln0  >  4- ^  the  asymptotic  expansion  of  the  first 
term  describes  the  head  waves.  The  condition  I'or  the  existence 
of  head  waves  is  that  one  or  more  of  the  -p;  be  less  than../^,. 

The  number  of  head  waves  which  exist  at  k  particular  offset  is 
determined  by  the  *P.  which  satisfy 
0-=^-P^-iSU1©^  '  jl. 

Associated  with  each  head  wave  there  is  a  critical  cone 

) 

a  critical  distance 
P  M 

and  a  travel  time 


_-£x  y 
^  vvr  Hpt-  z. 


V  _El  /  Q« 


The  asymptotic  expansion  of  U  can  be  expressed  In  the 


form 


J.2  i-”'  ^ 


(13) 


The  Tr  sre  normalized  travel  times  for  the  least-time  paths. 

J 

This  result  Indicates  that  the  response  function  and  all  Its 
derlvates  are  continuous  for  (j  =  The  first  term 

in  (13)  comes  from  the  saddle  points.  This  term  describes  the 
high  frequency  behavior  of  waves  which  travel  least -time  reflection 
paths.  %  Is  give  by 


M 


E,  A, 


1-1  HR(4rSUT%)''^ 

The  form  of  depends  on  whether  the  observation  point 
is  off  or  on  a  critical  cone.  At  a  point  which  is  not  on  a 


critical  cone 
oO 


/yi-l  ) 

The  leading  coefficients  are 


(1^) 


ij/  =  _  (Sine) 

vU  __i  I 


f  _ 


sLno 


(3501% 


sin^e  V  is, 


<3, 
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The  are  determined  by 

^  Z_  Hr  q'  /L  Hr  Qi^  ^  “’"zL  ^ 


A  = 


_f-EiSt(£±i^ 


7 


02. 


l=- 1  l- 

Each  of  the  coefficients  in  (14)  is 


where  Q-^  -sin^O 

continuous  as  0  approaches  zero.  This  is  somewhat  surprising 

for  in  the  derivation  of  (l4)  the  Hankel  functions  were  replaced 

by  their  asymptotic  expansions  for  large  argument  and  as  0 

approaches  zero,  the  argument  approaches  zero. 

As  the  point  of  observation  approaches  a  critical 

cone,  the  derivatives  of  -P  )  diverge.  This  follows  from  the 

fact  that  -p  ('^  )  is  a  function  of  the  square  roots, 

(J=l,  ,l).  Differentiation  introduces  terns  which  contain 

l/( -f as  a  multiplicative  coefficient.  These  quanti- 
J" 

ties  are  evaluated  at  ^  =sin0  .  As  the  point  of  observation 
approaches  the  critical  cone  (  0_),  sln^->-f  and  l/( -f 

d  J  vT 

diverges.  Consequently,  the  asymptotic  expansion  is  not  uni¬ 
form  in  sin@  . 

The  significance  of  becomes  apparent  when  we  note 

that  -f  (sin0)  is  a  function  of  the  ordlnarj'^  plane  wave  reflection 

and  transmission  coefficients  for  the  least-time  reflection  path. 

1/2 

The  factor  l/(  ^2^)'  corrects  the  plane  wave  amplitude  for 

geometric  spreading.  This  factor  can  be  u£-=?d  to  compute  the 
effect  of  geometric  spreading  for  any  3 east -time  reflection 
path  in  an  axially-symmet ric,  stratified  system.  Tiie  asymp¬ 
totic  expansion  in  (13)  is  derived  on  the  assumption  that  the 
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source  radiation  field  is  isotropic.  If  the  source  fiela  is 
axially  symmetric  but  not  isotropic,  contains  an  additional 
factor  which  specifies  the  dependemce  of  the  source  field  on 
the  polar  angle . 

At  a  point  which  lies  on  a  critical  cone 
oo 

~  (15) 

The  saddle  point,  sin9  ,  coincides  with  the  branch  point  -T- .  In 

U 

order  to  use  Laplace's  method  to  obtain  the  asymptotic  expansion 
or  (11),  -f  (£  )  must  be  analytic  at  the  saddle  point.  Analyticity 
is  achieved  by  introducing  the  transformation  A 


one-to-one  correspondence  between  points  in  the  and  ^  planes 
is  achieved  by  cutting  the  ,)|f-plane  along  the  imaginary  axis 


between  l-P-  and 

J  J 


The  “XcCJ)  can  be  expressed  in  terms 

A 


of  -p {>f)  and  its  derivatives  evaluated  at  the  origin  on  the  right 
side  of  the  cut.  The  leading  coefficients  are 


-y  (j.  _  - 


..  ■_  2'/'^r(3/4)  ,0^311/8/' 4 


TT 


1/2 


3/2. 


The  must  be  evaluated  at  sin  9=  -P  >  It  is  Important  to  note 

0 

that  on  a  critical  cone  p  (3lnQ  = -p  )='y(J)’  We  refer  to  the 

I  J~  I 

leading  term  as  the  geometric  term  and  to  the  difference  between  the 


exact  spectrum  and  the  geometric  term  as  the  perturbation  term.  On 
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a  critical  cone  the  asymptotic  expansion  is  obtained  in  powers 

\/A 

of  'y  ,  elsewhere  the  expansion  is  obtained  in  powers  of  'J  . 
Consequently  the  perturbation  spectrum  must  exhibit  a  sharp 
peaking  of  the  high  frequency  components  near  the  critical  cones. 

The  second  and  thlixi  terms  in  (13)  describe  the  high 
frequency  behavior  of  the  primary  head  wave  and  the  secondaiy 
head  waves  respectively.  Again  by  introducing  the  transforma¬ 
tion  )f=(  we  make 

'  sJ 


analytic  at  the  branch  point 

J 

are 


(16) 


The  asymptotic  expansions 


The  leading  coefficients  are 


(17) 
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_ io/p 

Because,  of  the  terins  in  (l~^p)  ,  the  ^(j)  diverge  as  the 

point  of  observation  approaches  the  critical  cone.  The  same 
behavior  is  exhibited  by  the  ^  (for  the  least -time  reflection 
path).  The  geometric  factor  in  is  identical  with  that 

obtained  by  Heelan  (1953)* 


D.  The  Singular  Behavior  of  the  Generalized  Response 
Functions  at  Times  Associated  with  Least -Time  Paths 

The  high-frequency  character  of  the  spectrum  is  com¬ 
pletely  determined  by  the  non-analytlc  behavior  of  the  response 
function  at  the  .  If  the  response  function  were  analytic, 
the  spectrum  would  exhibit  a  completely  different  hlgn-f requency 
character  (an  analytic  function  decays  more  rapidly  with  fre¬ 
quency  than  Cl)  ^  where  N  Is  any  positive  Integer  no  matter  how 
large).  The  *3^  are  the  only  non-analytlc  points,  otherwise  there 
would  be  additional  contributions  to  (13). 

To  discover  the  nature  of  the  singularity  at  'll,  we 
apply  the  Inverse  Fourlei'  transform  to  the  geometric  term  and 
to  the  leading  terms  In  the  asymptotic  expansion  of  the  pertur¬ 
bation  term.  Each  of  these  terms  can  be  associated  with  a  whole 
family  of  time  functions.  All  members  of  a  family  possess 
spectra  which  exhibit  Identical  high  frequency  characteristics. 
Accordingly,  we  attach  significance  only  to  those  characteristics 
of  the  time  function  which  are  exhibited  by  all  members  of  the 
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family.  One  such  characteristic  is  the  highest  order  singularity. 
A  discontinuity  in  the  function  itself  is  regarded  as  a  higher 
order  singularity  than  a  discontinuity  in  slope,  etc.  The  high¬ 
est  order  singularity  associated  with  each  member  of  a  particu¬ 
lar  family  will  be  identical.  For  example,  if  the  function  is 
continuous  but  its  first  derivative  diverges  logarithmically  at 
Xj  every  member  of  the  family  will  be  continuous  and  will  have 
a  derivative  which  diverges  logarltlimlcally  at^^. 

The  functions  to  which  we  apply  the  inverse  Fourier 
transform  are  of  the  form 


—  ^  L  LCr/l'Tl 


(18) 


If  ^  ^0  and  1>  this  function  is  not  Fourier  transformable. 
However,  there  is  a  whole  family  of  time  functions  which  in  the 
high  frequency  limit  exhibit  the  dependence  on '~i  given  in  (l8). 
To  find  a  member  of  the  family  we  modify  the  low  frequency 
behavior  to  make  the  result  Fourier  transformable.  We  make  the 


modification  in  such  a  way  that  the  phase  and  amplitude  spectra 
are  respectively  odd  and  even  functions  of'Y.  One  way  of  modi¬ 
fying  (l8)  is  to  replace  ( 1*7  by  (ilY  +  OL )^where  Ct  is  positive 
real.  It  is  Important  to  note  that  in  what  follows  none  of  the 
results  depends  on  Of. 
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The  time  function  corresponding  to  the  modified  ^^('7^) 


can  be  obtained  by  convolving  the  time  functions  corresponding  to 


LC'Tf)  —  A^e 
e  J 


:  ^r/iTl 


and 


Ffr) 


— 


In  the  process  of  convolution, 


-h 


^(T)cos^  I 


j 


produces  the  same  results  as  L  (  O' )  and  may  be  considered  equiva¬ 


lent  to  it .  ^  ( T )  is  the  Dirac  delta  function.  The  time  func¬ 

tion  which  corresponds  to  P  is 

P(T).0,  T-T^, 


Li 


«0 


2'"V'^r1 

e  convol 

I 


Applying  the  convolution  theorem  gives 


^  AjC0S^Pm,T-T^ 


(19) 


TT 


— oO 


The  singular  character  of  the  response  function  at  each  7^  can  be 
obtained  by  substituting  the  Individual  terms  from  the  asymptotic 
expansion  into  (l8)  and  (19) • 


V 
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The  geometric  term  for  the  least-time  reflection  path 


is 


The  strengths  of  the  delta  function  and  (0'-%)  terms  are  deter¬ 
mined  by  the  real  and  imaginary  parts  of  the  function  of  the 
plane  wave  coefficients  respectively.  Inside  and  on  the  first 
critical  cone  )  Is  pure  real  and  the  (^-'X)  * 

term  is  absent.  Except  for  the  fac..or  which  corrects  for  geo¬ 
metric  spreading,  R('X  )  similar  to  a  result  derived  by  Arons 
and  Yennie  (1950,^  by  considering  the  reflection  of  a  plane  wave 
at  a  liquid-solid  interface.  Equation  (20)  shows  that  their 
result  is  also  valid  for  a  multiply  reflected  wave  which  is 
initially  spherical  in  shape. 

The  geometric  terms  associated  with  the  head  v/aves 
produce  finite  discontinuities  at  the  X  given  by 

w 


This  is  a  real  quantity  because  pure  Imaglna  v. 

,3| -j- =0(J/1 ),  ^1  I  =1,  It  may  be  worth  pointing  out  that  (21) 
does  not  predict  how  the  head  wave  amplitude  decays  with  distance 
but  only  how  the  discontinuity  decays  with  distance.  Under 
certain  conditions  the  two  are  equivalent  but  not  always.  Con¬ 
sider  the  response  produced  by  a  source  function  ^ 

We  have  shown  that  convolution  of  a  function  which  varies  linearly 
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with  a  physically  realistic  source  function  (which  satisfies  (2)) 
yields  zero.  Consequently,  if  on  eacn  side  of ,  the  impulse 
response  varies  linearly  over  an  Interval  T  which  Is  equal  to 
the  duration  of  >1',  then  (21)  describes  how  the  head  wave  ampli¬ 
tude  decays  with  distance.  Certainly  for  a  source  function  of 
finite  duration  this  condition  cannot  be  satisfied  In  the 


vicinity  of  a  critical  cone  but  It  may  be  approached  as  the 
distance  from  the  critical  cone  Increases, 


Convolving  with  the  step  discontinuities 

associated  with  the  head  waves  Introduces  a  factor  In  the 


exoress,.  un  for  the  head  wave  amplitude 


2y  dividing  ouux'Cti 


I'j-  //  /  /\\'^ 

function  by  we  male  the  total  energy  radiated  Into  the 


system  independent  of  .  It  follows  that  the  head  wave  ampll- 
tude  Is  proportional  to  ^/(lyi^jeven  though  the  energy  Input  is 
independent  of  "tjj  .  As  the  dominant  period  Increases,  the  energy 


In  the  .head  wave  increases.  If  this  behavior  were  continued  to 


very  long  periods,  the  energy  In  the  h'ad  wave  would  exceed  the 
total  energy  radiated  Into  the  system.  We  conclude  that  a 

low-frequency  dependence  of  the  head  wave  amplitude 
is  Indicative  of  the  failure  of  the  high-frequency  theory. 

The  singular  behavic r  of  the  response  computed  from 
geometric  theory  Is  sketched  In  Figure  2. 


-  107 


Next  we  Investigate  the  singular  behavior  of  the  per¬ 
turbation  function  at  the  .  The  singularities  are  sketched  in 

Figure  3-  At  a  point  which  is  not  on  a  critical  cone,  the 
reflected  wave  is  associated  with 


Inside  the  first  critical  cone^  is  pure  real  and  the  log¬ 
arithmic  term  is  absent.  On  a  critical  cone  the  reflected  wave 
is  given  by 


Pa)' 


A.si.nC  C 


(23) 


RCt)- 


'rm)\ 


A.  Sin 


where 


li  ^Tr/8 


Ae"-e  1CI)= 


2^h^(3/4)  f  -fj- 

,  'J2.  !  -  -  oyi5  1  ^ 


IT' 


c= 


r  df 


3/2 


ai* 


i-.  fCf-O 


C  =  P  f- 

)  f  (f+l)^ 


P  indicates  that  the  Cauchy  principal  value  is  required.  On 
tie  first  critical  cone  (j=l),  iif  ^  is  pure  real.  Con- 

sequently,  the  perturbation  function  vanishes  for  'R  To  •  On 


any  critical  cone  except  the  first,  the  perturbation  function 


108  - 


diverges  from  both  the  left  and  right.  On  the  critical  cones 
the  divergence  is  not  symmetric  about  %  while  everywhere  else 


it  Is  symmetric. 

The  head  waves  are  associated  with  discontinuities 
In  the  slope  of  the  perturbation  function  which  are  given  by 


(24) 


By  convolving  a  physically  realistic  source  function 
with  each  of  the  geometric  singularities  (e.g.,  the  step,  delta 
function,  and  (T—To)  *  term),  we  obtain  three  primary  wavelets. 
The  high-frequency  response  of  a  layered  system  consists  of  a 
superposition  of  the  wavelets  for  all  least -time  paths.  By 
time  delaying  the  wavelets  and  multiplying  by  the  proper  ampli¬ 
tude  factors,  we  can  construct  the  high-frequency  response  from 
three  basic  waveforms. 


E .  The  One-Interface  System 

For  a  system  which  consists  of  two  half-spaces  per¬ 
fectly  coupled  together  aJong  a  plane,  several  Interesting 
characteristics  of  the  spectrum  can  be  Inferred  directly  from 
(9)  without  resrrting  to  approximate  theory.  The  geometric 

parameters  are  indicated  In  Figure  4.  Let  V  =^V  and  H  =R 

K  *  K 

(the  distance  along  the  least -time  path  for  the  reflected  com- 
presslonal  (P)  wave).  The  spectrum  of  the  iv-sponse  function 
produced  by  waves  which  arrive  at  the  receiver  In  the  comp'^es- 


slonal  mode  is 
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.0^ 


—  I  -•YiicosSU)^^- 0'^'^ 

(RV|)  U  =  ^7^ (-yf  su)epp)-f(f)e  ^  -o  ^ 


oif  ^  'y=c  CO  ^ 


T  ■ 


O 

ABOVE 


The  spectrum,  considered  as  a  function  of  Y  ,  has  a  shape  that 
depends  only  on  the  angle  of  incidence,  .  Changes  in  source 
and  receiver  positions  which  do  not  change  ^pp  have  no  effect 
on  the  spectral  shape.  The  spectrum  plotted  as  a  function  of  CO 
is  displaced  toward  lower  frequencies  as  R  increases  and  toward 
higher  frequencies  as  R  decreases.  In  general,  the  shear 
spectrum  is  more  complicated.  However,  if  the  source  and 
receiver  are  situated  in  the  same  pl'^ne  (Z=H),  both  the  com- 
pressional  and  shear  spectra  can  ’•'e  put  in  the  form 

(Hv;)u=iT'  f  i(2Yf+»ne„p)-fti')e  ^  '  "dc 

J 

0 

ABOVE 


0U-|-v 
Vi  ' 


where  ^=1  for  the  compressional  component  and  ^=0  for  the 
shear  component.  From  (25),  we  infer  that  the  shape  of  the 
response  function  (considered  as  a  function  of  "b/r^Vj)  )  depends 
only  on  ^pp •  This  means  that  the  amplitudes  o.  the  head, 
reflected  and  Interface  waves  bear  a  fixed  relationship  to  one 
another  which  is  independent  of  position  on  a  particular  ray. 
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If  we  impose  certain  criteria  to  determine  a  cutoff 
on  the  usefulness  of  the  high  frequency  theory,  that  cutoff  must 
satisfy  a  relation  of  the  form 

Increasing  H  reduces  and  extends  to  lower  values  the  fre¬ 

quency  ratige  within  which  the  use  of  high  frequency  theory  is 
Justified.  In  near  zone  work  whose  objective  is  to  delineate 
the  layering,  we  are  not  interested  in  frequencies  below  a 
certain  value.  By  using  (26)  we  can  find  a  function,  H^(  $pp)* 
which  corresponds  to  the  low-frequency  cutoff.  This  function 
defines  a  fictitious  surface  lying  above  the  interface.  Above 
this  surface  the  use  of  high-frequerwy  theory  is  Justified 
over  the  entire  frequency  band  of  interest --but  between  this 
surface  and  the  interface  it  is  not . 

Cagnlard's  method  yields  an  expression  for  the  exact 
impulse  response.  Inside  the  critical  angle  the  response  func¬ 
tions  for  the  compresslonal  and  shear  components  can  be  easily 
separated  into  a  geometric  term  (which  is  Just  the  delta  func¬ 
tion)  and  a  perturbation  term.  The  perturbation  tenn  is  plotted 
in  Figure  5  for  V'/V^a=0,5(  =30° )  and  in  Figure  6  for  Vj/V2  =1*5 

(no  critical  angle).  ^  is  the  Poisson’s  ratio.  In  each 
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figure  the  perturbation  function  is  plotted  for  several  values 
of  ^rp  ‘  The  compressional  functions  are  plotted  on  the  left 
side  under  P'  and  the  shear  functions  are  plotted  on  the  right 

s 

side  under  P  .  ^  =t/(H/V|)-  Except  where  indicated  otherwise, 

the  separation  between  divisions  on  the  abscissa  is  the  same 
for  all  graphs  on  a  particular  figure  and  is  indicated  in  the 
upper  left  graph.  The  long-time  divergent  behavior  is  charac¬ 
teristic  of  the  response  function  for  a  single  generalized 
transmission  path.  The  total  perturbation  response  is  obtained 
by  adding  the  compressional  and  shear  components  with  the 
appropriate  time  delay.  The  resultant  function  approaches  a 
constant  in  the  long-time  limit.  Convolution  with  a  source 
function  which  has  no  static  C''>mponent  produces  a  function  which 
tends  to  zero  in  the  long-time  limit. 

In  Figure  5  the  perturbation  functions  are  feature¬ 
less  except  for  the  discontinuity.  The  discontinuity  develops 
into  a  sharp  spike  as  approaches  the  critical  angle.  The 

asymptotic  theory  predicts  that  the  discontinuity  approaches 

infinity  as  0pp— >  0^  and  that  right  on  the  critical  cone 

'  ^  -3/-f 

the  initial  behavior  is  given  by  (T-Tpp) 

The  condition  for  the  existence  of  a  Stoneley  wave 
is  not  satisfied  in  Figure  6,  The  strong  phase  at  grazing 
incidence  is  a  pseudo-interface  wave.  Waves  of  this  type  have 
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been  studied  extensi  v?ely  by  St  rick  (1959)>  Phlnney  (1961)  and 
Gilbert  and  Laster  (1962).  There  are  no  terns  in  the  asymptotic 
expansions  for  either  high  frequency  or  large  horizontal  distance 
which  would  reveal  the  existence  of  pseudo -Interface  waves.  These 
waves  are  associated  with  zeros  of  the  denominator  In  the  gener¬ 
alized  coefficients  which  lie  on  a  lower  (non-permlsslble )  sheet 
of  the  Rlemann  surface.  The  pseudo-interface  wave  propagates 
with  a  velocity  which  exceeds  the  smallest  body  wave  velocity 
and  continually  radiates  energy  Into  one  or  both  of  the  adjacent 
media.  The  velocity  of  the  true  Interface  wave  is  always  less 
than  the  minimum  body  wave  velocity  and  there  Is  no  net  flux 
of  energy  normal  to  the  Interface.  In  Figure  6  the  velocity 
of  the  pseudo-interface  wave  is  Intermediate  between  the  two 
shear  velocities.  Therefore,  the  pseudo-interface  wave  produces 
motion  In  medium  one  but  there  is  no  net  flux  of  energy  normal 
to  the  Interface  which  is  directed  Into  medium  one. 

Amplitude  and  phase  spectra  have  been  computed  for 
the  perturbation  functions  of  Figure  5-  Individual  frequency 
components  )  are  plotted  against  the  angle  of  Inci¬ 

dence  In  Figures  7  and  8  for  the  ccmpresslonal  and  shear  com¬ 
ponents  respectively.  The  pronounced  minima  exhibited  by  the 
higher  frequency  components  occur  at  angles  where  the  discon¬ 
tinuity  goes  through  zero.  As  7^  decreases  the  minimum  becomes 
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less  pronounced  and  shifts  to  larger  angles.  The  lov;er  fre¬ 
quency  components  show  very  little  dependence  on  ^pp  •  higher 

frequencies  increase  in  amplitude  very  rapidly  as  0pp  approaches 
the  critical  angle.  Recall  that  asymptotic  theory  predicted  a 
sharp  peaking  of  the  high  frequency  components  on  each  critical 
cone . 

The  geometric  term  is  represented  by  a  delta  function. 
Therefore,  at  each  angle  of  incidence  all  frequencies  have  the 
same  amplitude.  The  amplitude  is  indicated  by  the  dashed  curves 
in  Figures  7  and  8.  Spectral  components  which  Intersect  a  ver¬ 
tical  line  above  the  dashed  curve  are  strongest  in  the  pertur¬ 
bation  term.  A  spectral  component  is  the  crossover  frequency 
at  that  angle  where  it  Intersects  the  dashed  curve.  The  cross¬ 
over  frequency  sets  a  lower  limit  on  the  frequency  range  over 
which  the  use  of  geometric  theory  is  Justified.  For  the  com- 
presslonal  component  the  crossover  frequency  lies  between  .25 
and  .50  over  nearly  the  entire  range  of  ©pp,  while  for  the 
shear  component  the  crossover  frequency  varies  from  about  1.0 
to  infinity.  The  reason  for  the  difference  is  that  the  geo¬ 
metric  term  for  shear  goes  to  zero  at  normal  incidence  and  Just 
Inside  the  critical  angle  while  the  geometric  term  for  the  com- 
pressional  component  does  not  go  to  zero.  At  angles  where  the 
geometric  term  vanishes,  it  is  the  discontinuity  in  t  ?  pertur¬ 
bation  function  which  determines  the  high  frequency  behavior  of 
the  total  response. 
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Let  T  be  the  duration  of  a  source  function  and 
its  dominant  frequency.  Whether  the  perturbation  term  can  be 
neglected  in  the  interval  T  following  the  onset  depends  on  the 
relation  between  the  crossover  frequency  and  .  When  is 
large  compared  to  the  crossover  frequency,  the  perturbation 
term  can  be  neglected.  Outside  the  Interval  T  following  the 
onset  the  situation  is  more  complicated.  The  perturbation 
term  may  contain  additional  phases  (e.g.,  pseudo -interface  waves) 
which  do  not  overlap  with  the  reflected  {diase.  Although  these 
phases  are  diminished  in  amplitude  by  increasing  4^  ,  they 
probably  cannot  be  neglected  unless  they  Interfere  with  larger 
amplitude  geometric  phases  reflected  from  other  Interfaces  or 


fall  outside  the  dynamic  range  of  the  recording  system. 
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VI.  ANALYSIS  OF  A  PROGRAM  FOR 
COMPUTING  THEORETICAL  SEISMOGRAMS 
FOR  A  MULTILAYERED  MEDIA 

The  "matrix  method"  (Haskell,  1953)  Is  in  principle  a 
method  for  the  numerical  computation  of  theoretical  seismograms. 
The  method  is  analyzed  with  the  objective  of  detemining  the 
effect  of  source  depth  on  the  response  of  an  elastic  layered 
medium  in  the  time  interval  preceding  the  arrival  of  the  unat- 
tenuated  normal  modes.  The  mathematical  model  is  based  on  an 
idealized  physical  model  of  a  vacuum  half-space  overlying  a  set 
of  liquid  and  solid  elastic  layers,  with  all  liquid  layers  above 
the  solid  ones.  The  layers  are  homogeneous,  isotropic,  and 
perfectly  coupled  together  at  the  interfaces.  The  point  source 
and  the  point  receiver  may  be  located  in  any  layer. 

The  mathematical  technique  is  the  common  one  of 
assuming  that  the  stresses  and  displacements  may  be  derived 
from  a  set  of  potential  functions  (one  shear  potential  and  one 
compressional  potential  for  each  layer),  expressing  these  poten¬ 
tials  in  terms  of  their  Fourier  transforms  and  then  using  the 
matrix  method  to  determine  the  transform  of  the  response,  i.e., 
the  transforms  of  the  stresses  and  displacements.  In  theory, 
the  actual  physical  stresses  and  displacements  may  be  found  by 
inverting  these  transfonns  by  double  integration  over  the  first 
quadrant  of  the  plane  of  k  (wave  number)  and  w  (circular  fre¬ 
quency),  the  transform  parameters.  However,  in  general  there 
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is  E  finite  positive  number  such  that  the  transforms  are 

singular  for  some  curves  in  the  region  where  w  /K  <  and 
regular  for  w/K>Cg^  .  is  the  maximum  phase  velocity  for 

unattenuated  normal  mode  propagation.  In  this  work  the  inte¬ 
gration  is  constrained  to  the  regular  region. 

Other  investigators  have  used  a  similar  technique  but 
have  confined  their  attention  to  the  location  of  the  singular 
curves  (dispersion  curves)  and  the  contribution  to  the  total 
r*esponse  of  these  singular  curves,  i.e.,  they  have  investigated 
the  normal  modes,  which  yield  the  most  significant  contributions 
to  the  response  at  large  epicentral  distances.  Additionally, 
there  are  certain  technical  differences  in  the  definitions  in 
the  matrix  method  itself,  resulting  in  the  fact  that  in  our 
formulation  the  matrix  which  relates  the  stress-displacement 
vector  at  the  bottom  of  a  layer  to  the  one  at  the  top  of  the 
layer  is  real  for  all  non-negative  K  and  w  .  The  point  com- 
pressional  source  is  also  introduced  here  by  real  quantities, 
free  of  singularities. 

The  formulas  were  programmed  in  Fortran  for  the 
7094,  with  the  exception  that  the  complex  subroutine  was  not 
used.  Computation  of  the  transformed  response  vector  V  for 
each  K  and  w  in  the  case  of  a  vacuum  over  a  solid  half-space 
requires  twenty  milliseconds.  Each  additional  solid  layer 
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requires  about  five  milliseconds .  For  each  w,  v  is  computed  for 
a  set  of  equally  spaced  K  values  running  from  zero  up  to  the 
cut-off  valve  of  w'/C^  .  The  w  values  are  also  equally  spaced, 

starting  with  w  about  ten. 

The  transformed  response  V  for  a  vacuum  over  a  solid 
half -space  exhibits  sharp  peaking  at  the  point  K  =  (  <x  is 

the  compresslonal  velocity).  In  an  Interval  of  the  K  values  of 
10“3,  one  component  rises  from  zero  to  6500.  Since  these  values 
are  to  be  Integrated,  very  fine  sampling  of  the  v  functions 
is  necessary  in  the  K  direction.  Because  it  is  desirable  that 
the  source  pulse  be  sharp,  a  wide  range  of  w  values  is  also 
required,  necessitating  sampling  the  v  functions  at  over  a 
million  points  in  the  (  K,w) -plane.  This  sampling  alone  would 
require  five  and  one-half  hours  on  the  computer,  not  including 
the  double  Integration  wiilch  must  fellow. 

Certain  difficulties  occur  even  during  the  computation 
of  V  .  Elements  of  the  matrices  Involve  quantities  such  as 
exp  )y  where  h  is  the  thickness  of  a  layer  or 

combination  of  layers.  This  number  is  frequently  outside  the 
range  of  the  computer  (with  the  standard  floating  point  subroutine) 
for  many,  indeed,  most,  problems  of  seismologlcal  Interest.  Also, 
expressions  such,  as 

"I  K^-CW/p<)^  "I  (v//|$)^  h 

appear,  which  lead  to  large  round-eff  errors  for  large  h  or  K  . 
These  problems  have  been  experienced  by  others,  who  have  circumvented 
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them  by  using  the  very  low  w  values  and  the  small  k  values 
which  are  more  applicable  in  the  normal  mode  theory.  Theee 
p-wblems  of  number  size  and  lound  off  can,  of  course,  be  solved 
by  1. -creasing  the  word  length  of  the  computer,  either  mechani¬ 
cally  or  by  special  programming. 

The  effect  of  restricting  the  w  and  K  values  to  the 
region  of  regularity  vv/k>C^  was  partially  investigated  by 
exam.inlng  the  effect  of  such  a  truncation  on  the  potential  func¬ 
tion  for  a  point  compressional  source  in  an  infinite  homogeneous 
isotropic  medium.  It  was  discovered  that  the  distortion  caused 
in  this  potential  by  eliminating  the  region  w/K<Cg  was  not 
great  if  the  receiver  was  on  the  same  horizontal  as  the  source, 
but  was  much  greater  if  the  receiver  was  on  the  same  vertical 
as  the  source.  Because  of  this  distortion,  it  would  be  diffi¬ 
cult  to  determine  in  the  seismogram  what  role  was  played  by  the 
layering  in  the  system  and  what  role  was  played  by  the  placing 
of  the  receiver. 

It  was  concluded  that  though  the  matrix  method  is 
conceptually  a  simple  numerical  method  for  solving  this  complex 
problem,  in  practice  too  much  computation  is  required  to  realize 
even  the  result  distorted  by  truncation,  w.hich  would  be  diffi¬ 
cult  to  Interpret.  Certainly,  the  expense  is  currently  too 
great  to  justify  use  o"'  this  method. 
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The  matrix  method  has  some  value  as  a  theoretical  tool 
here,  particularly  regarding  location  of  the  curves  of  singulari¬ 
ties.  The  following  conclusions  were  reached  by  a  simple  exami¬ 
nation  of  the  appropriate  mathematical  expressions; 

1)  The  only  branch  points  of  v  in  the  case  of  the 
point  compressional  source  are  those  involving  the  compressional 
and  shear  velocities  in  the  last  layer  (the  lower  half -space), 
and  this  is  true  for  any  intermixed  system  of  liquid  and  solid 
layers.  Jardetzky  (1953)  proved  this  for  solid  layers. 

2)  If  all  layers  are  liquid  with  compressional 
velocities  ,  then  the  phase  velocity  C  of  any  unattenuated 

Q 

normal  mode  satisfies  min  <x.  ^  where  is  the  com- 

pressional  velocity  of  the  half -space.  In  particular,  if 
min  ^  i  there  are  no  unattenuated  normal  modes. 

3)  In  the  case  of  liquid  layers  over  a  solid  half- 

space,  if  min  ,  where  tx*  is  the  Rayleigh 

velocity  of  the  half -space,  there  is  at  least  one  unattenuated 
normal  mode  with  a  phase  velocity  smaller  than  . 
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VII.  THE  REFRACTED  ARRIVAL  PROM  A  LAYER 

High-frequency  geometric  ray  theory  is  used  to  investi¬ 
gate  the  refracced  arrival  from  a  high-speed  layer  embedded  in 
an  infinite  medium.  The  geometric  parameters  are  indicated  in 
Figure  1.  The  top  of  the  layer  is  at  a  depth  H  below  the  source 
and  at  a  depth  Z  below  the  receiver.  The  Z  axis  passes  through 
the  source  and  is  normal  to  the  Interface.  The  plane  Z=0  coin¬ 
cides  with  the  top  of  the  layer.  The  radiation  field  is  axially 
symmetric  about  the  Z  axis.  E  is  the  layer  thickness  and  p  is 
the  range.  V  is  the  compressional  velocity,  nr  is  the  shear 
velocity,  and  b  is  the  density.  Subscript  1  refers  to  the 
infinite  medium  and  subscript  2  refers  to  the  layer. 

Refraction  along  a  layer  is  considerably  more  compli¬ 
cated  than  refraction  along  an  interface  between  two  semi¬ 
infinite  media.  To  demonstrate  this,  let  T  be  the  duration  of 
the  refracted  arrival  from  a  semi -infinite  medium  and  let  "tTfs  be 
the  travel  time.  Note  that  tu  is  independent  of  layer  thickness. 
The  results  for  a  semi -infinite  medium  are  not  valid  if  wa/es 
which  are  multiply  reflected  within  the  layer  or  other  wave  types 
arrive  in  the  interval  ^  "t  ^  .  Multiple  reflections 

distort  the  refracted  arrival  when  T  is  large,  the  layer  is  thin, 
or  the  offset  is  large.  When  any  of  these  conditions  is  satis¬ 
fied,  the  individual  f^ases  are  not  completely  resolved  and  the 
character  and  amplitude  of  the  refracted  arrival  are  determined 
by  the  way  in  which  the  individual  phases  interfere  with  one 


another. 
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When  we  attempt  to  take  Into  account  all  the  waves 
which  are  multiply  reflected  within  the  layer,  we  immediately 
encounter  tne  problem  of  degeneracy.  The  simplest  case  of 
degeneracy  arises  when  the  ray  crosses  the  layer  once  in  the 
compressional  (P)  mode  and  once  in  the  shear  (S)  mode.  In 
Figure  1  we  see  that  chere  are  two  such  rays  which  arriv^'  at 
exactly  the  same  time  and  are  therefore  degenerate.  Thi  se  two 
rays  differ  in  the  sequence  in  which  the  P  and  S  legs  of  the 
path  are  traversed.  This  type  of  degeneracy  is  entirely  differ¬ 
ent  from  the  accidental  degeneracy  v/hlch  can  occur  when  the 
/elocities  and  the  layer  thickness  are  suitably  related.  It  is 
a  consequence  of  the  fact  that  two  different  types  of  propaga¬ 
tion  occur  in  an  elastic  medium.  Degeneracy  occurs  only  for  rays 
which  cross  the  layer  in  both  modes. 

To  reduce  this  idea  to  quantitative  terms,  we  intro¬ 
duce  the  following  notation: 

h  =  number  of  crossings  in  the  P  mode, 
yn  -  number  of  crossings  in  the  S  mode. 

Specifying  n  and  Tf]  does  not  determine  a  particular  ray  but 
rather  a  family  of  degenerate  rays.  The  number  of  rays  in  the 
family  (i.e,  the  degeneracy)  is 


(n-t-fn) ! 
n!m ! 


(1) 
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Even  for  small  n  and  Tf\  ,  the  degeneracy  can  be  quite  large.  If 
highly  degenerate  events  are  truly  significant,  a  way  must  be 
found  for  removing  the  degeneracy- -otherwise  the  ray  theory 
becomes  extremely  inefficient  from  a  computational  standpoint. 
Also,  because  all  the  degenerate  waves  arrive  at  the  same  time, 
only  the  composite  event  has  physical  significance. 

For  the  problem  of  an  embedded  layer,  a  simple  method 
has  been  found  for  removing  degeneracy.  This  method  permits  us 
to  express  the  total  amplitude  of  a  degenerate  event  (Aj^^  )  in 
terms  of  a  simple  series  which  contains  either  rf)  or  H  +1  terms 
(whichever  is  smaller).  Effectively  we  replace  the  computation 
of  the  amplitudes  of  individual  rays  by  the  evaluation  of  nf} 

or  f) +1  terms  in  a  series.  Furthermore,  the  use  of  the  series 
avoids  the  complicated  bookkeeping  required  to  Insure  that  all 
members  of  the  degenerate  family  have  been  included. 

We  study  the  refracted  arrival  produced  by  a  source 
which  radiates  a  spherically  symmetric  com;-re3Slonal  wave.  The 
vertical  component  of  the  particle  velocity  in  the  direct  wave 
is 


is  the  distance  between  the  source  and  receiver.  The  source 
function,  )^,  is  time-limited  and  merely  expands  or  contracts 
with  ‘tjji  but  does  not  change  shape.  The  amplitude  spectrum  of 
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is  assumed  to  peak  at  a  dominant  frequency  =  1/^^  .  At  dis¬ 

tances  which  are  large  compared  with  the  dominant  wavelength, 
the  direct  wave  propagates  without  change  of  shape  and  the  shape 
is  determined  by  the  source  function  ^  is  the  radius  of 

the  assumed  spherical  cavity  at  the  source  and  is  the  peak 
pressure  applied  to  the  cavity  wall.  The  factor 
is  introduced  to  make  the  total  energy  input  independent  of  . 

Beyond  a  certain  minimum  range  which  depends  on  the 
time  interval  of  interest,  the  only  reflected  waves  which  can 
contribute  to  the  refracted  arrival  are  those  whose  phase  veloci¬ 
ties  approach  the  compressional  velocity  in  the  layer  at  large 
distances.  According  to  geometric  ray  theory,  the  only  reflected 
waves  which  satisfy  this  condition  are  those  which  cross  the 
layer  at  least  once  in  the  compressional  mode  (D^l).  The  high- 
frequency  asymptotic  representation  for  the  composite  reflection 
when  \fl  ^  1  is 


n-sl  • 


Snm  A 


is  just  the  sum  of  products  of  the  plane  wave  reflection 


and  transmission  coefficients  and  is  the  normalized  travel 


time.  The  time  variation  in  each  reflected  wave  which  contributes 


to  the  refracted  arrival  is  determined  by  the  source  function  )f . 
determines  the  reduction  in  amplitude  produced  by  geometric 
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spreading  and  is 


^nn.~  (■ 


i+# 


»/  / 


))  ^ 


)  w 


where  is  the  angle  of  Incidence  at  the  upper  interface  of  rays 
which  cross  the  layer  V)  times  in  the  compressional  mode  and  ttl  times 
in  the  shear  mode.  For  1,  =sin'“%^g^. 

The  reflected  waves  that  do  not  cross  the  layer  in  the 
compressional  mode  (H  =0  )  generate  head  waves.  The  head  waves 
that  propagate  with  a  phase  velocity  contribute  to  the  refracted 
arrival.  The  direct  P  wave  generates  the  ordinary  head  wave  indi¬ 
cated  by  the  segment  ED  in  the  upper  diagram  in  Figure  2.  It  also 
generates  a  transmitted  compressional  wave  (P,  P^  )  and  a  transmitted 
shear  wave  The  transmitted  wavefronts  intersect  at  the 

interface  as  long  as  the  phase  velocity  in  the  direct  wave  exceeds 
.  VRien  the  phase  velocity  in  the  direct  wave  drops  below  ,  the 
wavefronts  separate.  In  this  process  an  internal  head  wave  is 
generated  (BD).  This  wave  travels  in  the  direction  normal  to  BD 
with  the  velocity'll. 

The  lower  diagram  demonstrates  what  happens  when  the 
transmitted  shear  wave  and  the  head  wave  reflect  off  the  base  of 


the  layer.  BQ  is  the  reflected  part  of  the  head  wave.  The  inci¬ 
dent  shear  wave  (PC)  generates  a  reflected  compressional  wave  (QF) 
and  a  reflected  shear  wave  (PBA).  These  two  reflected  wavefronts 
begin  to  separate  at  a  point  on  the  lower  interface  where  the  phase 


velocity  in  the  incident  shear  wave  drops  below  .  A  new  head 
wave  is  generated  whose  wavefront  coincides  with  BQ.  Each  time 
the  unconverted  shear  wave  (n=o)  is  reflected  back  into  the  layer, 
a  new  head  wave  is  generated  whose  wavefront  coincides  with  the 
wavefront  of  head  waves  generated  by  earlier  shear  reflections. 

The  high-frequency  asymptotic  representation  for  head 
waves  which  propagate  with  a  phase  velocity  is 


UJ  I  pj  HY, ’  H  ’  H  'JT/  I H/Vii 


where 

The  time  variation  is  determined  by  the  integral  of  the  source 
function.  At  large  offsets  the  amplitude  decays  like  ( P/H)  , 
and  the  relative  amplitudes  are  determined  by  the//^  .  The/Z^yj 
depend  on  the  compresslonal  velocity  ratio,  the  density  ratio, 
and  the  Poisson’s  ratios.  The  wavefront  of  each  head  wave  is 
tangent  to  the  reflected  wave  which  generated  it  at  the  critical 


distance  p 


m 


'Tyyi  is  the  normalized  travel  time.  The  head  wave 


phases  are  separated  by  equal  time  Intervals 
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which  do  not  depend  on  offset.  This  fact  suggests  that  if  the 
refracted  arrival  consisted  of  head  waves  alone,  the  condition 
"tj  would  lead  to  constructive  interference  when  the  all 

have  the  same  polarity. 

Comparison  of  (3)  and  (5)  shows  that  in  the  high  fre¬ 
quency  limit  the  reflected  waves  predominate  over  the  head  waves. 
The  complete  expressjon  for  the  refracted  arrival  is 

C -"Z 

I  -•  rn=*o^2)’*’ 

is  the  travel  time  of  the  earliest  arriving  reflected  wave 

which  is  incident  beyond  the  critical  angle  The  minimum 

range  at  which  (7)  can  be  applied  is  determined  by 


AT, 


where  A  T  is  the  time  Interval  of  interest.  Equation  (7)  indi¬ 
cates  that  except  for  the  factor  (  b|V|  /\  ),  the  response 

expressed  as  a  function  of  i;/(H/V,  )  depends  on  the  three  geometric 
parameters  p/H,  E/H  and  Z/H,  and  on  the  dominant  wavelength-to- 
thickness  ratio  through  the  parameter 

ii, 

m  H  V  E  A  H  ;  ■ 

For  high  frequencies  the  multiply -reflected  waves  deter¬ 
mine  the  characteristics  of  the  refracted  arrival.  The  travel - 
time  curves  for  reflected  waves  for  which  ffl  =0,1,2  are  plotted  in 
Figure  3,  The  difference  in  travel-time  between  the  reflected 
wave  and  the  refracted  arrival  is  plotted  along  the  ordinate.  The 
range  is  plotted  on  a  logarithmic  scale  to  emphasize  the  fact  that 
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all  the  travel -time  curves  for  a  particular  HI  value  approach  the 
same  asymptote  at  large  ranges  regardless  of  the  f)  values .  The 
travel -time  difference  between  waves  which  have  different  T(\  values 
approaches  a  finite  value  at  large  offsets  which  is  a  multiple  of 
is  defined  in  (6),  where  it  is  shown  that  2T^  is  the 
time  interval  between  head  wave  arrivals . 

This  clustering  of  the  travel -time  curves  leads  us  to 
refer  to  waves  with  the  same  rn  value  as  a  particular  order.  If 
the  significant  amplitudes  in  each  order  are  confined  to  a  limited 
range  of  H  values,  and  if  T  (source  duration)  is  less  than  , 
the  refracted  arrival  at  large  ranges  consists  of  a  sequence  of 
events.  Each  event  is  associated  with  a  particular  fT)  value  and 
propagates  with  a  phase  velocity  which  is  very  near  the  compresslonal 
velocity  in  the  layer.  The  amplitude  of  each  event  is  obtained  by 
summing  the  amplitudes  of  all  rays  which  have  the  same  in  v ilue  as 
follows ; 


m 


A  study  of  many  cases  reveals  that  for  p/H^20  and  E/H-^.l, 


(8) 


At  large  ranges  the  head  wave  contributions  decay  like  (p/H)~  . 
This  means  that  the  reflected  waves  determine  the  character  of 
the  refracted  arrival  at  large  distances. 
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The  'Yjyj  determine  the  relative  amplitudes  of  the  events 
associated  with  the  different  nO  values.  The  depend  on  the 
compresslonal  velocity  ratio,  the  density  ratio,  and  the  Poisson's 
ratios.  A  study  of  many  cases  reveals  that  the  higher  orders  (ni^'l) 
have  amplitudes  which  are  at  least  ten  per  cent  of  the  amplitude 
for  rn  -C  when  either  the  compresslonal  velocity  contrast  is 
large,  the  Poisson's  ratio  in  the  layer  is  near  0.25,  or  the 
Poisson's  ratio  in  the  infinite  medium  approaches  Q 5* 

Equation  (8)  cannot  be  used  when  the  individual  wavelets 
(corresponding  to  the  dlfierent  fl  values)  are  wholly  or  partially 
resolved.  The  degree  of  resolution  Increases  as  the  offset 
decreases,  the  layer  gets  thicker,  or  the  dominant  frequency 
increases.  In  Figure  4,  the  amplitude  of  each  wave  is  plotted 
along  the  ordinate  and  the  difference  between  the  arrival  time 
of  the  reflected  wave  and  the  onset  of  the  refracted  arrival 
along  the  abscissa.  Here  we  consider  only  multiples  which  do 
not  cross  the  layer  in  the  shear  mode  (ni-O).  Each  curve  is 
drawn  for  a  particular  offset.  The  leftm.ost  dot  on  each  curve 
gives  the  amplitude  and  time  difference  for  the  wave  which  crosses 
the  layer  twice  in  the  compresslonal  mode  (n=2).  As  we  move 
along  each  curve  from  left  to  right,  R  Increases  in  increments 


of  two. 
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The  Infomatlon  displayed  In  Figure  4  predicts  the 
phenomenon  of  shingling.  This  effect  Is  generally  explained  In 
terms  of  normal  mode  propagation  In  which  the  phase  velocity  exceeds 
the  group  velocity.  As  the  range  Increases,  peaks  and  troughs 
move  forward  through  the  envelope  which  defines  the  refracted 
arrival.  In  this  process,  the  amplitude  of  the  first  extremum 
decreases  and  is  eventually  lost  in  the  noise.  At  this  offset, 
a  later  extremum  Is  selected  to  define  the  time-distance  curve. 

At  each  offset  where  an  extremum  Is  lost,  there  Is  a  discontinuity 
In  the  time-distance  curve  and  a  new  shingle  Is  added  correspond¬ 
ing  to  a  later,  larger  amplitude  extremum.  In  Figure  4,  we  note 
that  as  the  range  Increases,  the  time  delay  decreases  for  each 
n  value.  The  wavelet  associated  with  each  reflection  moves  for¬ 
ward  In  time  with  respect  to  the  onset  and  the  amplitudes  of  the 
earliest  arrivals  (D -2  and  4)  decrease.  As  the  range  increases 
the  amplitudes  of  the  later  arrivals  Increase  (the  open  dot 
indicates  this  effect  for  H  =8),  and  the  number  of  arrivals  within 
a  fixed  time  interval  Increases.  When  the  wavelets  are  partially 
resolved,  this  causes  additional  wiggles  to  appear  on  the  tall 
of  the  first  event.  The  m '.vement  of  each  wavelet  forward  toward 
the  onset,  the  decrease  In  amplitude  of  the  early  extrema 
(associated  with  small  If]  values)  and  the  appearance  and  increase 
In  amplitude  of  the  later  extrema  (associated  with  large  fl  values) 
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act  to  create  the  shingling  effect.  Clearly,  in  shingling  the 
dui’atlun  of  the  individual  wavelets  must  not  be  so  short  that 
they  appear  as  discrete  events.  On  the  other  hand,  the  duration 
must  not  be  so  long  that  all  resolution  is  lost. 

The  reflected  waves  and  head  waves  superpose  to  foimi 
the  refracted  arrival.  The  actual  response  function  is  obtained 
by  (a)  convolving  the  source  function  with  the  sequence  of 
impulses  which  represent  the  reflected  waves  and  (b)  convolving 
the  integral  of  the  source  function  with  the  sequence  of  impulses 
which  represent  the  head  waves  (as  indicated  in  Equations  (3), 

(5)  and  (7)).  The  choice  of  the  source  function  is  not  completely 
arbitrary.  It  must  be  chosen  so  that  if  the  medium  which  con¬ 
tains  the  source  were  infinite,  each  point  would  return  to  its 
original  state  a  finite  time  after  the  arrival  of  the  direcc  wave. 
This  requires  that  the  static  components  of  the  stress  and  strain 
be  finite.  The  Laplace  transform  of  the  radial  stress  in  a 
spherically  symmetric  compressional  wave  js 


In  order  for  the  static  component  to  be  finite 

iimCtr )(  —  A j  ^—2  . 

S-^o 

As  a  consequence  of  (9),  must  have  at  least  two  axis  crossings. 
We  also  require  that  the  particle  velocity  and  acceleration  be 
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continuous.  The  function 

satisfies  all  these  requirements.  The  first  factor  determines 
the  envelope.  The  duration  of  the  function  is  T-Ki::^  .  The  condi¬ 
tion  2  forces  the  low-frequency  behavior  to  satisfy  (9).  The 
source  functln  and  its  Integral  are  plotted  in  Figure  5  K=4. 

The  source  function  is  symmetric  about  "t =2.  The  integral  of 
the  source  function  is  antisymmetric  about  “t /tj  =2.  Under  cer¬ 
tain  conditions  we  can  dlstlnguis’.i  between  reflected  and  head 
wave  contributions  to  the  refracted  arrival  on  the  basis  of  sym¬ 
metry. 

Figures  6-b  show  how  the  character  of  the  refracted 
arrival  varies  with  the  dominant  frequency  at  three  offsets .  In 
each  figure  the  dimensions  are  fixed.  One  division  along  the 
abscissa  is  equal  to  the  dominant  period.  As  we  move  down,  the 
dominant  period  Increases,  the  dominant  wavelength  increases, 
and  decreases.  Because  the  dominant  period  changes,  the 

time  scales  cannot  be  compared  directly. 

The  refracted  arrival  for  an  infinitely  thick  layer 
would  have  the  antisymmetric  fom  of  the  head  wave  and  teminate 
at  the  fourth  division  on  all  traces.  Figure  6  shows  three  dis¬ 
tinct  events  on  the  traces  designated  E/y^^  =7j5j  and  4.  These 
events  correspond  to  in  =0,1,  and  2.  As  E/^j  decreases,  the 
later  events  move  toward  the  first  event.  Interfere  with  it,  and 


become  lost  in  it . 
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Changes  in  E/  produce  significant  changes  in  charac¬ 
ter  for  p/H  values  of  3-5  and  10,  but  relatively  minor  changes  for 
p/H=100.  At  large  ranges  the  dominant  frequency  is  not  high  enough 
to  partially  resolve  the  individual  waves  which  contribute  to  the 
first  event.  At  small  ranges,  the  individual  waves  are  separated 
by  larger  time  intervals  and  the  change  in  character  is  a  con¬ 
sequence  of  interference. 

Increasing  the  dominant  period  causes  the  individual  waves 
to  overlap  to  such  an  extent  that  only  a  single  event  is  discern¬ 
ible  and  reduces  the  amplitude  of  the  reflected  viaves  relative  to 
the  head  waves.  This  explains  why  the  last  trace  on  each  figure 
approaches  the  antisymmetric  shape  associated  with  the  head  wave. 

Each  trace  is  normalized  so  that  the  maximum  amplitude 
plots  at  the  same  value.  The  true  maximum  amplitude  in  the  refrac¬ 
ted  wave  train  is  plotted  in  Figure  9  as  a  function  of  the  dominant - 
wavelength-to-thickness  ratio.  Each  curve  is  drawn  for  a  particular 
range.  As  the  range  increases,  a  relative  maximum  develops  which 
migrates  toward  the  vertical  line  marked  have  shown  that 

at  large  ranges  the  waves  arrive  in  groups  which  are  associated 
with  different  fn  values  and  are  separated  by  the  time  interval 

.  "tr^  is  the  condition  for  constructive  interference  between 

groups.  The  weak  minimum  and  maximum  at  p  /H=100  are  also  a  con¬ 
sequence  of  interference  between  the  groups. 
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As  the  dominant  frequency  Increases,  the  individual 
waves  in  the  group  begin  to  separate  from  one  another,  and  the 
amplitude  varies  in  a  complex  way.  Destructive  interference 
between  individual  waves  produces  the  deep  minimum.  At  sufficiently 
high  frequencies,  the  individual  waves  are  completely  resolved, 
and  each  curve  must  approach  an  asymptote  with  slope-1/2.  At  very 
low  frequencies  the  head  wave  is  dominant,  and  each  curve  must 
approach  an  asymptote  with  slope  1/2.  The  increase  in  amplitude 
at  low  frequencies  is  indicative  of  the  failure  of  the  high- 
frequency  theory.  If  this  increase  were  continued  to  very  low 
frequency,  the  energy  in  the  refracted  arrival  would  exceed  the 
total  energy  radiated  into  the  system.  We  do  not  attach  sig¬ 
nificance  to  dominant -wavelength-to -thickness  ratios  greater 
than  unity.  Consequently  our  analysis  is  restricted  to  thick, 
high  velocity  zones . 

Figures  10-12  Illustrate  the  effect  of  range  for 
E/ ==10,  7  and  4.  On  each  figure,  E/^is  fixed  and  the  dominant 
period  does  not  change  from  trace  to  trace.  Each  trace  is  nor¬ 
malized  so  that  the  maximum  amplitude  plots  at  the  same  value. 

The  nonna3ized  range  (  p/H)  changes  in  increments  of  0.5.  The 
group  of  lines  that  start  in  the  upper  left  corner  of  Figure  10 
shows  how  individual  peaks  and  troughs  move  forward  with  respect 
to  the  onset  and  decrease  in  amplitude  as  the  range  increases.  The 
decrease  in  amplitude  would  be  considerably  more  pronounced  if  true 


amplitudes  had  been  plotted.  The  second  group  of  lines  shows  how 
additional  peaks  and  troughs  appear  on  the  tall  of  the  fir^t  event 
as  the  offset  increases.  These  phases  are  associated  with  the 
later  arriving,  higher  H  values  which  increase  in  amplitude  with 
range.  The  third  group  of  lines  shows  that  the  second  event  (inf\=l) 
also  exhibits  shingling. 

In  Figure  11,  the  phase  velocity  is  slightly  greater  for 
the  third  peak  than  it  is  for  the  earlier  extrema.  Shingling  is 
exhibited  only  over  a  limited  range.  The  second  trough  and  third 
peak  on  the  first  trace  can  be  correlated  across  the  entire  record, 
but  the  phase  velocity  difference  disappears  for  P/I1>“7*  In 
Figure  12,  there  is  no  shingling.  These  observations  indicate 
that  shingling  is  range-limited.  The  maximum  range  decreases  as 
the  dominant  period  increases,  and  for  sufficiently  long  dominant 
periods,  shingling  does  not  occur.  When  shingling  is  observed, 
it  indicates  that  the  high  velocity  zone  is  thick.  Increasing 
the  layer  thickness,  decreasing  the  dominant  period,  and  decreasing 
the  layer  depth  all  act  to  increase  the  effective  time  interval 
between  wavelets  (l.e.,  the  degree  of  resolution)  and  thereby 
extend  the  maximum  range  over  which  shingling  can  occur.  Range - 
limited  shingling  is  distinctly  different  from  the  range -independent 
shingling  associated  with  normal  mode  propagation.  The  Importance 
of  shingling  lies  in  its  potential  to  determine  layer  ti. 'ckness. 
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VIII,  REFLECTION  AND  TRANSMISSION  OF 
PLANE  COMPRSSSIONAL  WAVES 

The  reflection  and  transmission  of  a  plane,  mono¬ 
chromatic,  compressional  (P)  wave  at  a  plane  interface  between 
two  semi-infinite  media  is  studied  in  great  detail  in  Vela 
Scientific  Report  No.  4  (AFCRL  64-205).  That  report  contains 
extensive  tables  of  the  energy  flux  ratio  in  the  vertical  direc¬ 
tion  and  of  the  relative  phase  of  the  vertical  displacement  for 
each  of  the  reflected  and  transmitted  waves.  Our  purpose  in 
treating  this  problem  is  three -fold: 

1.  To  provide  the  scientific  community  with  a  compre¬ 
hensive  set  of  data  which  demonstrates  the  effect  of  systematically 
varying  the  compressional  velocity  ratio,  the  density  ratio,  and 
the  Poisson's  ratios  over  a  wide  enough  range  to  include  cases 

of  interest  in  geophysical  applications,  seismic  modeling,  and 
velocity  measurements. 

2.  To  arrive  at  certain  general  conclusions  about  the 
effect  of  the  parameters  on  the  partition  of  the  incident  flux  among 
the  different  waves  and  on  the  angular  dependence  of  the  flux  in 
each  wave . 

3.  To  clarify  the  physical  significance  of  Knott’s 
equation  when  the  angle  of  incidence  exceeds  the  critical  angle. 
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The  four  parameters  which  we  select  to  describe  the 
media  are  the  compressional  velocity  ratio  V21,  the  density  ratio 
R021,  and  the  Poisson's  ratios  SIGl  and  SIG2.  The  incident  P 
wave  approaches  the  Interface  through  medium  1  and  is  transmitted 
into  medium  2.  The  notation  V21  is  used  to  represent  the  ratio 
of  the  compressional  velocity  in  medium  2  to  the  compressional 
velocity  in  medium  1.  Similarly  R021  represents  the  density  ratio. 

The  following  range  of  parameters  has  been  investigated; 


V21- 

-0. 

.25, 

0.50, 

0, 

.75, 

1, 

o 

o 

1, 

.50, 

ro 

o 

3.0  and  4.0j 

R021- 

-0. 

.33, 

0.50, 

0, 

,80, 

0, 

.90, 

1, 

1.10, 

,  1.20,  1.50, 

2. 

,0  and  3.0 

• 

> 

SIGl- 

-0. 

,10, 

0.20, 

0, 

.25, 

0. 

.30, 

0, 

,40 

and  0, 

.50; 

SIG2- 

-0, 

,10, 

0.20, 

0, 

.25, 

c, 

.30, 

0, 

,40 

and  0, 

,50. 

A  material  with  a  Poisson’s  ratio  of  0.50  is  a  fluid. 

The  flux  ratios  and  phase  angles  have  been  tabulated  for  all  com¬ 
binations  of  the  above  parameters  except  where  the  fluid  velocity 
is  higher  than  the  P  wave  velocity  of  the  adjacent  solid  or  where 
the  fluid  density  is  greater  than  that  of  the  solid.  In  addition, 
we  have  studied  the  cases  where  V21=0.10,  R021=0.0005,  3IG2=0.5 
and  SIGl  varies  from  0.10  to  0,50  as  given  above.  These  cases 
describe  a  wave  which  is  Incident  upon  an  air  boundary  from  a 
solid  or  liquid.  And  conversely,  we  have  treated  the  air  wave 
which  is  incident  upon  a  solid  or  liquid  boundary.  Here  V21=10, 
R021=2000,  SIG1=0.5  and  3IG2  varies  over  the  range  given  above. 
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For  each  of  some  two  thousand  different  combinations  of 
parameters,  the  energy  flux  ratios  and  the  relative  phases  of  the 
vertical  displacement  have  been  computed  at  five  degree  Intervals 
for  angles  of  incidence  ranging  from  0  to  85  degrees.  When  a 
critical  angle  is  present,  these  quantities  are  also  computed  at 
the  critical  angle  and  at  one  degree  intervals  within  five  degrees 
of  the  critical  angle. 

The  computations  Indicate  that  when  one  medium  is  a 
fluid  (or  air),  a  large  amount  of  shear  (S)  energy  is  generated 
in  the  solid  medium.  In  many  cases  this  conversion  has  an 
efficiency  of  better  than  50  per  cent  over  some  tens  of  degrees 
of  the  angle  of  incidence. 

Analytical  expressions  are  presented  for  the  reflection 
and  transmission  coefficients  for  the  vertical  component  of  the 
particle  displacement.  Each  coefficient  can  be  put  in  the  fonn 

i.  0{&)  I 

Q=A(0  )0  ,  where  &  is  the  angle  of  Incidence.  (j)  {  &  ) 

determines  the  phase  shift  between  the  incident  and  reflected 
(or  transmitted)  waves  and  is  listed  in  the  tables.  When  Q  is 
less  than  the  critical  angle,  all  the  coefficients  are  real 
(^0  or  Tf  ),  and  Q  determines  the  ratio  of  the  vertical  displace¬ 
ment  in  the  reflected  (or  transmitted)  wave  to  the  vertical  dis¬ 
placement  in  the  incident  wave.  When  Q  is  greater  than  the 
critical  angle,  all  the  coefficients  are  in  general  complex,  the 
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incident  and  reflected  (and  transmitted)  waves  are  not  in  phase, 
and  Q  no  longer  measures  the  ratio  of  the  vertical  displacements 
(for  the  ratio  is  time -dependent ) . 

Some  interesting  relations  are  obtained  by  evaluating 
the  analytical  expressions  for  the  coefficients  at  the  critical 
angles.  At  the  critical  angle  for  the  transmitted  P  wave 
(  =  sln'v, /V2_),  the  vertical  displacement  in  the  transmitted 

P  wave  vanishes  and  the  condition  for  perfect  reflection  with  no 
conversion  and  no  transmission  is 


,  QJi ,  and  p.^  are  respectively  the  compressional  velocity, 

’  >  ■th 

Shear  velocity,  and  density  in  the  i  medium.  When  this  condi¬ 
tion  is  satisfied,  all  the  energy  in  the  incident  P  wave  goes 
into  the  reflected  P  wave  and  the  vertical  displacement  undergoes 
a  l8C®  phase  change.  At  the  critical  angle  for  the  transmitted 
S  wave  (  =sin~'  V, /O^),  the  condition  for  perfect  reflection  is 
satisfied  if  either 

(2) 

or 

I  /  Fx 

)■ 


(3) 
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Each  of  these  conditions  Insures  that  no  energy  goes  into  the 
reflected  S  wave.  At  and  beyond  the  critical  angle  for  S,  the 
net  flux  (the  instantaneous  flux  integrated  over  a  period)  in 
each  of  the  transmitted  waves  vanishes.  There  may  be  an  instanta¬ 
neous  flux  in  the  transmitted  P  and/or  S  wave,  but  when  Integrated 
over  a  period,  it  averages  to  zero.  Therefore,  (2)  or  (3)  insures 
that  the  reflected  P  wave  is  the  only  wave  which  carries  a  net  flux 
of  energy  away  from  the  interface.  On  the  other  hand,  (l)  insures 
perfect  reflection  of  both  the  instantaneous  and  net  fluxes.  This 
distinction  between  the  instantaneous  and  net  fluxes  is  basic  in 
the  discussion  of  the  requirements  for  continuity  of  the  normal 
component  of  flux  across  the  interface. 

The  instantaneous  value  of  the  normal  component  of 
flux  in  each  wave  can  be  expressed  in  terms  of  the  A  and  ^ 
associated  with  that  wave.  Inside  the  critical  angle  {  ) ,  the 

vertical  and  horizontal  components  of  displacement  in  each  wave 
are  either  in  phase  or  18O®  out  of  phase,  the  particle  motion  is 
linear,  and  the  direction  (but  not  the  magnitude)  of  the  energy 
flux  is  independent  of  the  time.  If  the  "transmitted"  P  wave  is 
studied  for  9  >■  or  the  "transmitted"  S  wave  is  studied  for 

we  find  that  the  vertical  and  horizontal  displacements 
are  out  of  phase  by  90'’,  the  particle  motion  is  elliptical  and 
the  direction  (and  magnitude)  of  the  energy  flux  oscillates  with 
time.  Actually,  at  each  point  the  flux  direction  reverses  every 
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quarter  cycle,  and  at  each  Instant  the  flux  direction  reverses 
every  quarter  wavelength  along  the  interface.  Therefore,  beyond 
the  critical  angle  for  a  transmitted  wave,  the  net  flux  in  that 
wave  vanishes  and  the  energy  is  cycled  back  and  forth  between  the 
two  media. 

Inside  the  critical  angle,  ,  the  conditions  for  the 
CO!  tinuity  of  the  normal  component  of  the  Instantaneous  flux  and 
of  the  normal  component  of  the  net  flux  are  identical  and  are 
expressed  by  Knott's  equation.  This  equation  contains  the 
associated  with  each  wave.  Beyond  the  critical  angle,  Knott's  equa¬ 
tion,  by  itself,  is  not  sufficient  to  Insure  continuity  of  the 
instantaneous  f lux--actually  three  equations  are  required.  The 
first  is  obtained  from  Knott's  equation  by  setting  the  A  for  the 
transmitted  P  wave  equal  to  zero  when  (or  the  A's  for 

both  the  transmitted  P  and  S  waves  equal  to  zero  when  ). 

This  modified  Knott's  equation  guarantees  the  continuity  of  the 
net  flux.  The  other  two  equations  are  just  the  real  and  Imaginary 
parts  of  the  complex  equation  obtained  by  substituting  AG  for 
A  in  Knott's  equation.  Inside  the  critical  angle  the  terms  in 
Knott's  equation  may  be  interpreted  as  either  Instantaneous  or  net 
flux  ratios.  Beyond  the  critical  angle  the  terms  in  Knott's 
modified  equation  can  be  Interpreted  only  as  net  flux  ratios. 
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IX.  SOME  EXPERIMENTS  ON 
INVERSE  FILTERING  OF  SEISMIC  RECORDS 

Our  original  research  proposal  tentatively  suggested 
an  investigation  of  the  application  of  inverse  filtering  to  the 
determination  of  seismic  energy  source  depth.  This  was  intended 
to  be  an  offshoot  of  our  own  geophysical  exploration  research 
program,  which  has  for  several  years  Included  a  study  of  general 
filtering  problems  together  with  the  limitations  imposed  by  back¬ 
ground  noise.  Inverse  filtering,  which  has  also  been  called 
inverse  convolution  or  deconvolution,  is  in  fact  very  intimately 
linked  with  concepts  of  optimum  filtering,  such  as  those  intro¬ 
duced  by  Norbert  Wiener  (19^9).  In  recent  years  there  have  been 
many  communications  pertinent  to  the  subject.  Including  the 
tutorial  article  by  G.  L.  Turin  (i960)  and  those  by  R.  B.  Rice 
(1962),  G.  Kunetz  (1961),  M.  M.  Backus  (l96i),  J.  d'Hoerane  (1962), 
M.  R.  Foster,  et  al  (1962),  and  £.  A.  Robinson  (1957). 

The  possible  application  of  inverse  filtering  to  the 
detection  and  identification  problem  is  perhaps  best  visualized 
with  the  aid  of  the  schematic  diagram  in  Figure  1 .  It  is  quite 
similar  to  the  general  problem  posed  by  G.  L.  Turin  (196O)  in 
his  ^‘igure  3.  In  the  problem  of  Interest  here  it  is  considered 
that  the  design  of  the  processing  filter  is  such  that,  when  com¬ 
bined  with  the  conventional  receiver  filter,  it  forms  a  composite 
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optimum  inverse  filter.  This  composite  filter  should  be  aimed 
at  yielding  an  output  v;hlch  is  a  best  possible  representation 
of  the  seismic  sour’ce  conditions.  It  should  ideally  be  a  true 
inverse  filter  for  the  transmission  characteristic  of  the  earth, 
but  in  practice  j.t  must  be  optimized  in  some  sense  to  take  into 
account  the  estimated  signal-to-noise  ratio.  The  concept  of 
such  an  optimum  inverse  filter  is  fairly  straightforward,  but 
in  order  to  design  it  well,  it  is  evld.  n.tly  necessary  to  know 
the  relevant  transmission  characteristic  of  the  earth  in  con  • 
slderable  detail.  In  fact,  it  is  lack  of  detailed  knowledge  of 
this  transmission  characteristic  that  mainly  hinders  the  applica¬ 
tion  of  inverse  filtering  to  the  detection  and  identificavluij 
problem.  Some  earth  transmission  characteristics  which  have 
become  reasonably  well  known  in  rejent  years  are  the  low  fre¬ 
quency  phase  dispersions  of  the  first  modes  of  the  Love  and 
R'-’lelgh  surface  waves.  At  large  ranges  these  dispersi  iS  are 
very  considerable.  Hence  it  is  not  surprising  to  finu  that 
inverse  filtering  has  already  been  applied  to  such  surface  waves 
by  K.  Akl  (i960),  J.  N.  Brune,  et  al  (i960),  and  Y.  Sato  (1955> 
1956).  Further  research  into  this  application  of  inverse  filter¬ 
ing  to  surface  wave  data  is  currently  being  directed  by 
Professors  M.  Ewing  and  F.  Press  under  Project  VELA  contracts. 

We  have  not  attempted  to  perform  such  investigations  at  Calresearch 
as  there  would  seem  to  be  little  point  in  mainly  duj  Heating  the 
efforts  at  the  Lament  Geological  Observatory  and  the  California 
Institute  of  Technology. 
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Some  detailed  knowledge  of  long  range,  or  third  zone, 
direct  body  wave  transmission  through  the  earth  has  recently  been 
obtained  by  the  United  Kingdom  group  using  their  Pole  Mountain 
array.  The  Russian  nuclear  event  of  February  2,  1962,  was 
recorded  with  a  reasonable  signal -to -noise  ratio  ^.n  the  0.7  to 
1.7  cps  band.  The  main  significant  arrival  was  a  "lonesome  P" 
pulse.  A  reproduction  of  this  direct  P  pulse,  as  recorded  by 
one  of  the  "quietest"  Pole  Mountain  seismometers,  is  shown  in 
Figure  2A.  The  pulse  is  extremely  simple  in  character.  In  fact, 
if  it  is  assumed  t.hat  the  time  behavior  of  the  source  function 
is  fairly  accurately  known,  th'^^n  it  may  be  considered  that  the 
recording  in  Figure  2a  effectively  specifies  the  composite  filter 
of  the  Willmore  seismometer  and  the  89*’  path  between  the  location 
of  the  Russian  bomb  and  Pole  Mountain.  On  this  basis,  it  might 
at  first  sight  seem  reasonable  to  design  an  inverse  to  this 
composite  filter.  This  inverse  filter  could  then  be  applied  to 
some  other  direct  P  event,  following  a  similar  path  of  about 
39®,  in  order  to  yield  its  approximate  source  function.  There 
is  at  least  one  stiong  argument  against  such  a  course  of  action. 
It  is  evident  from  the  simple  pulse  shape  of  Figure  2a  that  the 
composite  transmission  path  and  seismometer  filter  is  essentially 
band-pass  in  nature  with  relatively  slight  phase  distortion.  It 
is  this  simplicity,  of  course,  which  makes  third  zone  detection 
schemes  look  so  attractive.  An  optimum  inverse  filter  might  well 


164 


nari’ow  the  direct  P  pulse  in  Figure  2a,  but  practical  signal-to- 
noise  limitations  will  normally  prevent  a  dramatic  shortening 
and  simplification  of  the  pulse.  Nevertheless,  the  importance 
of  the  identification  problem  is  such  that  even  a  marginal 
improvement  in  the  over-all  filtering  system  may  be  worth  striv¬ 
ing  for. 

Following  reasoning  somewhat  analogous  to  this,  but 
oriented  more  towards  emphasizing  first  oreaks,  J.  W.  Birtill 
of  the  United  Kingdom  Atomic  Energy  Authority  conducted  some 
relevant  analog  Inverse  filtering  experiments  (private  communi¬ 
cation).  In  order  to  supplement  J.  W.  Birtill ’s  analog  experi¬ 
ments,  some  digital  studies  were  later  instigated  by  us,  using 
a  slightly  modified  version  of  a  Calresearch  7090  program.  It 
must  be  emphasized  that  these  studies  were  far  from  exhaustive, 
the  lack  of  detailed  data  on  slgnal-to-nolse  ratios  of  the 
records  being  a  major  restriction.  Further  research  on  this 
topic  should  definitely  include  a  more  detailed  knowledge  of 
realistic  signal -to-noise  ratios  and  would  possibly  be  pursued 
more  rapidly  by  personnel  with  more  immediate  access  to  contem¬ 
porary  data  from  seismic  arrays . 

In  our  simple  investigation  we  assumed  that  the  effec¬ 
tive  source  function  for  the  Russian  bomb  was  doublet  in  displace¬ 
ment  (see  Figure  2b).  The  basis  for  this  was  our  belief  that. 
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at  a  large  range  from  a  detonation  within  an  infinite  medium,  the 
displacement  would  be  approximately  Impulsive.  We  considered  that 
the  presence  of  the  free  surface  Immediately  above  the  detonation 
would  generate  a  subsequent  impulse  of  the  same  strength  but  with 
reversed  polarity.  If  the  actual  surface  stresses  lie  in  a 
grossly  non-linear  region,  this  concept  of  ci  doublet  source  func¬ 
tion  may  be  substantially  incorrect  and  this  possibility  should 
not  be  forgotten. 

The  approximate  power  spectra  of  the  hypothesizea 
doublet  source  and  the  Russian  bomb  event  of  Figure  2  are  shown 
in  Figures  3  and  4  respectively.  The  ratio  of  their  amplitude 
spectra  is  shown  in  Figure  5  and  the  difference  between  their 
phase  spectra  is  given  in  Figure  6. 

By  definition  the  amplitude  and  phase  characteristics 
in  Figures  5  and  6  specify  the  inverse  filter  which  is  required 
to  change  the  Russian  bomb  event  (Figure  2a)  into  the  idealized 
doublet  source  function  (Figure  2b).  In  the  construction  of 
this  inverse  filter,  however,  it  is  only  sensible  to  note  the 
following;  although  Figure  4  specifies  a  spectrum  for  the 
Russian  bomb  event  between  0.0  and  3-0  aps,  it  is  evident  that, 
in  all  probability,  there  is  a  poor  signal -to-noise  ratio  outside 
the  regj-on  of,  say,  0.3  to  2.0  cps.*  This  means  that  the  filter 


*Much  of  the  noise  is  undoubtedly  digitizing  noise.  Unfortunately, 
the  ratio  of  signal-to-ambient  ground  motion  has  not  been  reliably 


estimated . 
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of  Figure  5^  which  has  a  very  strong  transmission  of  high  frequen¬ 
cies,  will  mainly  pass  noise  when  applied  to  the  Russian  event. 

It  would  be  illogical  to  utilize  such  an  inverse  filter  and,  as 
we  noted  earlier,  the  estimated  signal-to-noise  ratio  should  be 
taken  into  account  in  optimizing  this  filter.  However,  because 
we  had  no  reliable  estimate  of  the  signal-to-noise  ratio,  it  was 
decided  to  modify  the  spectrum  of  Figure  5  by  merely  truncating 
it  outside  some  pass  band  whose  high-  and  low-cut  frequencies  can 
be  specified  arbitrarily.  Approximate  impulse  responses  of  some 
resulting  inverse  filters  are  shown  in  Figures  7-10,  together 
with  their  truncation  frequencies. 

The  effect  of  applying  these  inverse  filters  to  the 
Russian  bomb  event  is  shown  in  Figure  11.  Attention  should  be 
given  only  to  those  pulses  synchronized  in  time  with  the  original 
Russian  event.*  It  is  fairly  evident  that  the  main  effect  of 
narrowing  the  pass  band  of  the  inverse  filter  is  to  increase  the 
effective  width  of  the  main  output  pulse.  The  output  pulse  is, 
in  each  case,  assyrnetrlcal  and  is  essentially  a  band  limited 
approximation  to  the  hypothesized  doublet  source . 

*The  surrounding  ripple  is  due  to  the  extremely  sharp  high- 
frequency  cut  off  of  each  inverse  filter  The  outer  pulses  are 
caused  by  undersampling  the  associated  Inverse  filter  character¬ 
istics  with  the  digitizing  rate  of  only  0.1  cps .  This  admittedly 
detracts  from  the  study  but  it  was  not  considered  worthwhile  to 
reprocess  the  data  wit.h  a  more  rapid  frequency  sampling  rate. 
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The  result  of  applying  the  inverse  filters  of  the 
Russian  bomb  to  the  French  bomb  of  May  1,  1962,  (fortunately  also 
at  an  epicentral  range  of  89°)  is  shovm  in  Figure  12.  It  is 
perhaps  significant  to  note  that  the  signal-to-ncis^^  ratio  was 
poorer  for  the  French  bomb.  This  may  account  for  the  rather 
large  oscillations  of  about  2.0  cps  introduced  with  the  0.3  to 
2.0  cps  Inverse  filter.  In  any  case  it  should  be  noted  that  our 
practical  inverse  filters  do  not  cause  gross  instabilities  when 
applied  to  the  French  bomb  record. 

As  a  final  step  the  inverse  filters  were  applied  to  an 
earthquake  from  the  Ionian  sea  (April  10,  I962).  The  results  are 
reproduced  in  Figure  I3  and  do  not  exhibit  any  particularly 
significant  features .  What  originally  looked  like  a  seismogram 
still  continued  to  look  like  a  seismogram  after  the  inverse 
filtering  operations. 

In  conclusions,  then,  we  have  reported  on  the  very 
simple  inverse  filtering  experiments  that  we  have  applied  to 
bomb  and  earthquake  data.  The  results  are  very  much  as  expected 
and  are  neither  dramatic  nor  dismaying.  We  do  not  plan  to  con¬ 
tinue  with  this  phase  of  research.  It  is  to  .;e  emphasized  that 
any  future  research  on  this  topic  should  incorporate  reliable 
estimations  of  signal-to-noise  ratios  and  would  possibly  best 
be  performed  by  organizations  with  more  immediate  access  to  con¬ 
temporary  data  from  sei.-mic  arrays. 


Figure  1.  Schematic  diagram  of  main  factors  in  the  detection  problem. 
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FIGURE  7. 
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